Radiative Hydrodynamics

Here is how the internal energy in the grid can change due to radiation:

\frac{\partial (\rho e)}{\partial t} = - \kappa_R \rho (4\pi B-cE)

where B=B(T)=\sigma T^4 is the bolometric Planck function for a blackbody (BB), and 4 \pi B=a T^4 is the energy output due to BB radiation. cE is the radiative energy density in the grid. This equation shows that when 4 \pi B > cE,

\frac{\partial \rho e}{\partial t}<0

which is interpreted as the matter losing energy via BB radiation. That is, the internal energy of that zone will decrease, having been transferred into radiative energy. This equation also tells us that when 4 \pi B<cE, there is more energy in the radiation field than in the BB, and so it gets absorbed by the matter. This causes the internal energy to increase,

\frac{\partial \rho e}{\partial t}>0

Changes to the internal energy lead to changes in temperature, and thus the next time step the amount of radiative energy from the BB will change (recall B=B(T)).

Now, since the matter and the radiation are coupled in this way, the equation that governs the radiative energy in the grid is the inverse of the internal energy. Additionally, the radiative energy can diffuse through the grid, so there is an extra term for diffusion:

\frac{\partial E}{\partial t}= \nabla \cdot (\frac{c \lambda}{\kappa_R \rho}) \nabla E + \kappa_p \rho(4\pi B-cE)

Note the coupling term comes in with a '+' sign now (rightmost term on the RHS), and the diffusion term (left term). Note also that how strongly the matter and radiation couple depends on the opacity.

Uniform ambient evolution

The rad. energy in the grid at t=0, assuming no 'sources' is given by the temperature field through the term:

E=4\pi B(T) = aT^4

In the code the constant out in front is given by "scalerad". Checking a uniform ambient medium with this expression for Erad, shows that Erad at t=0 is given by scalerad*T⁴. Thus, any gas with finite temperature is producing radiative energy in the grid through black body radiation. Since

\frac{\partial E}{\partial t} \propto \nabla \cdot \nabla E + (4 \pi B - cE)