The lack of interpolation points at a free surface and the corresponding increase in interpolation errors in an SPH calculation, dictates the addition of a boundary condition when a radiative surface condition is needed. This poses a number of problems, the form of the boundary condition and how to define the free surface.
For astrophysical problems the boundary can be defined approximately as the surface where radiation escapes freely. At the surface of a star the mean free path of the radiation becomes large because the opacity becomes small and therefore the diffusion approximation fails. In practice stellar structure calculations preserve the fiction that the diffusion approximation applies because the bulk of the star is not sensitive to the details of the treatment of the surface. To treat the surface in our SPH calculations we need to add a term to equation (11) so that particles near the free surface radiate freely.
Let us assume we have a particle radiating freely in space and it is cooling
like a black body sphere of radius 
. The rate at which it cools
per unit volume is
which can be written as
where 
is the opacity and we have assumed that 
.
The mean free path for a gas is given by 
, which in our case
is 
. Therefore if the particle radiates as a black body, the rate
at which it cools is
This is not a surprising results, and is obvious from a dimensional analysis of the divergence of the flux. This is the form that our boundary term will take so that equation (1), without a source term will become
where 
is a dimensionless scalar function that is zero
more than 
interior to the surface. If we integrate equation
(18) over the entire volume we get
Since 
is only nonzero within 
of the surface we can expand
Q and T in a taylor series around 
, the vector to the surface.
This is defined so that 
where
is the outward normal to the surface.
Expanding around the surface we get
where 
and we have assumed 
.
If we assume q is independent of where we are on the surface then equation
(20) becomes
where the integral is now over the surface and 
.
Therefore, equation (21) produces the correct radiative boundary condition if
Figure 3: A vector plot in the X-Y plane of 
as defined by equation (23) for a radiating equilibrium
figure with variable mass particles
The surface boundary function 
is easily constructed
by first defining the vector 
by,
or an alternate definition is given by
In both cases the integrand is an odd function of 
. If the
material within 2h of 
is evenly distributed then the
integral will be zero. If the material within 2h of 
is
not evenly distributed the integral will be nonzero reaching a maximum
at a free surface. Figure (3) is a plot of 
(defined by equation (23)) on all the particles within h
of the X-Y plane for a radiating equilibrium figure using variable mass
particles.
For equilibrium figures equation (23) or (24) will give a good approximation of the normal to the surface. Unfortunately for domains containing density discontinuities the functions above will be nonzero at the discontinuities making the definition of the free surface problematical.
Figure 4: Comparison of the analytic solution of the density with the
SPH solution for a radiating star.
Figure 5: Comparison of the analytic solution of the temperature with the
SPH solution for a radiating star.
Figures (4) and (5) show the density and temperature
profiles respectively for an SPH model of a radiating star using 3544
particles. The energy generation term was 
,
the equation of state of a perfect gas was assumed with 
and a Kramer's opacity law 
was used. The
choice for the energy generation was based on the desire for a star
radiative throughout, and without an extended envelope (see [10]).
As constant mass particles where used an extended
envelope would have regions with densities too low to be modelled.
From figure (4) we can see that even with the modest energy
generation used only one half of the radius can be modelled.
The radiating star was constructed using the full set of hydrodynamic equations with the addition of a damping term. Over time the particles relaxed to the low energy equilibrium figure shown in figures (4) and (5).
The outer envelope of low density material can be modelled by increasing the number of particles used, thus lowering the minimum density at each point, or by using variable mass particles.
