With the inclusion of a radiative heat flux term the energy equation can be written in the form
where S is the entropy/unit mass, T the temperature, 
the
radiative flux and 
is a possible energy source or sink term.
In the diffusion approximation the flux can be written as
where Q is a function of the opacity, temperature and density.
The standard SPH approximation for hydrodynamic equations is to replace spatial derivatives of functions with derivatives of the interpolation kernel W. An example of this would be the following symmetric SPH approximation of equation (2) (see [8], for a more complete description of standard SPH approximations),
where 
is the mass associated with the interpolation point j and
.
The interpolated value of the flux vector
from equation (3) can then be used in the SPH approximation of
equation (1) to get
This is essentially the method adopted by Lucy [5] in his study of binary fission. Computationally this is a straight forward adaption of the SPH solution of first derivatives, applied twice to calculate second derivatives. Unfortunately this method has a number of serious drawbacks. It is extremely sensitive to the interpolation point distribution, and radiating surfaces are not calculated correctly.
As a test of the method, Lucy constructed a three dimensional equilibrium model of a polytrope. He found that the calculated divergence of the flux had errors an order of magnitude greater than the pressure force calculation. This magnification of the error is due to the double interpolation from the particles. Any error in the flux calculation, due to particle disorder, or too few particles is compounded in the calculation of the divergence of the flux.
It has been argued that the SPH formulation implicitly incorporates an
adequate radiative boundary condition [5]. This boundary
condition is 
, 
as 
at
the free surface. For most astrophysical applications the free
surface boundary condition 
is adequate, for
radiative boundary conditions this is no longer true. In practice
because the number of particles contributing to the
interpolation is decreasing. As the number of points decreases the
solution is dominated by the shape of the interpolation kernel, not
the interpolated function. In the case of a radiating free surface it
is important that the radiating boundary condition be accurately
determined or applied. At a free surface this is impossible for SPH to
achieve (see [3], for a discussion of boundary conditions in
SPH). Without the incorporation of an explicit surface condition it
is impossible to construct a self gravitating equilibrium figure with
energy generation and heat diffusion [4]. The SPH
approximation is poor at the free surface and any reliance on an
implicit boundary condition is doomed to failure.
We present an alternate derivation of the divergence of the flux, one that does not incorporate a boundary condition. We then present a one dimensional test of the method demonstrating its robustness to particle dissorder. We then discuss the addition to the SPH approximation of equation (1) of a free surface radiating boundary condition and use the boundary condition to construct a three dimensional radiating equilibrium figure. Finally we present a stability analysis of the new SPH diffusion approximation.
