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),

cE=4\pi B ~(at ~t=0), and

\nabla E =0, we have

\frac{\partial E}{\partial t}=0

There is no change in the radiative energy over time. Now, other dynamics in the simulation, e.g. gravity, could change the temperature distribution in the grid and thus B(T). Only after this would happen, would we begin to see changes in E. Thus, while 4 \pi B = cE, the system is in radiative equilibrium (right?).

Jeans unstable gas

Imagine now starting with a uniform, Jeans unstable gas mass. Initially it is in radiative equilibrium (as in previous section), but within a freefall time, the gas will begin to collapse — becoming denser and hotter as it does (recall equation for internal energy has a gravitational energy term). This will lead to regions in the grid where 4piB>cE. This will increase the radiative energy field,

\frac{\partial E}{\partial t} \propto (4 \pi B - cE) > 0

and thus act to 'cool' the gas:

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

Now in the next timestep this radiative energy could either diffuse away, or stick around, depending on how optically thick the gas is (controlled by \kappa_R).

If the diffusion term is larger than the coupling term (which acts to increase Erad over the course of the infall by the conversion of gravitational energy into heat):

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

then, E will decrease faster than it is increasing (confused by the sign in the diffusion term), which implies:

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

i.e. the gas continues cooling. However, eventually the material should become dense enough that it becomes optically thick to the radiation, and the gas should begin to heat up. In order to get,

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

cE needs to get larger than 4piB. This can be achieved by the combined effect of increased BB radiation (through the compressional heating), in addition to slower diffusion. Thus, for this problem it seems you would want \kappa_R= \kappa_R (\rho). That is, in the early stages of collapse (low rho), the increased heat due to infall should be cooling through radiative losses. I.e., the collapse should remain isothermal. However, after a certain density is reached, the collapse should become adiabatic. This seems to be controllable through the diffusion term. I can't think of physical reasons why you might want to change \kappa_p. In what situations would you want more or less coupling?

Radiation with a source

Now we add a radiating source to the grid. Here is the initial condition:

Recall,

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

Thus, a zone adjacent to the source will acquire some E due to the sink (through the diffusion term). This is because of the gradient in E. There will be no contribution from BB radiation, as out there, E=B. Note, the strength of the diffusion depends on \kappa_R as well as the gradient in E. Here is the grid after the radiative time step:

The change in E is dominated by the diffusion term closer-in to the sink (strong gradient there), and further away it is dominated by BB radiation.

Now, depending on position, the internal energy will either decrease or increase. Recall,

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

So, close-in to the source, E>B, but further away (where diffusion was weaker), B>E. Thus, after a hydro timestep we have:

This means, depending on the details of the rad-hydro, there will be some regions of radiative heating and some of cooling. For this picture here (which is very rough btw — a diffusion wave, and coupling 'wave' may have different concavity than shown here), we have:

This process will continue until equilibrium has been reached — that is, no more gradients in E (and thus no more diffusion), and the energy in the radiative field matches the energy being produced by the blackbody (i.e. E=B). Recall, the ambient will not change because out there the only term in E is the BB contribution. That is, E=B out in the ambient, and thus, no changes occur there.