Physics of Radiation Transfer

tau = l \kappa=\frac{l}{\lambda_p
beta = \frac{u}{c
tau << 1	streaming limit
tau >> 1 \mbox{, } \beta \tau <<	static diffusion limit
tau >> 1 \mbox{, } \beta \tau >>	dynamic diffusion limit

Equations of Radiation Hydrodynamics

frac{\partial \rho}{\partial t} + \nabla \cdot \left ( \rho \mathbf{v} \right ) =

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P+\mathbf{G

frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = c G^

frac{\partial E}{\partial t} + \nabla \cdot \mathbf{F} = -c G_0

frac{1}{c^2} \frac{\partial \mathbf{F}}{\partial t} + \nabla \cdot \mathbf{P} = -\mathbf{G

where

E=\int_0^\infty d \nu \int d \Omega I(\mathbf{n}, \nu

mathbf{F}=\int_0^\infty d \nu \int d \Omega \mathbf{n} I(\mathbf{n}, \nu

\mathbf{P}=\int_0^\infty d \nu \int d \Omega \mathbf{nn} I(\mathbf{n}, \nu

and

G^0 = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right

\mathbf{G} = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right ] \mathbf{n

Simplifying assumptions

• If the flux spectrum of the radiation is direction-independent then we can write the radiation four-force density in terms of the moments of the radiation field

^0=\gamma [ \gamma^2 \kappa_{0E} + (1-\gamma^2)\kappa_{0F}]E-\gamma \kappa_{0P}a_RT_0^4-\gamma(\mathbf{v} \cdot \mathbf{F}/c^2)[\kappa_{0F}-2\gamma^2(\kappa_{0F}-\kappa_{0E})] - \gamma^3(\kappa_{0F}-\kappa+{0E})(\mathbf{vv}) : \mathbf{P}/c^

mathbf{G} = \gamma \kappa_{0F}(\mathbf{F}/c)-\gamma \kappa_{0P}a_RT_0^4 (\mathbf{v}/c) - \left [ \gamma^3( \kappa_{0F} - \kappa_{0E})( \mathbf{v}/c)E + \gamma \kappa_{0F}(\mathbf{v}/c) \cdot \mathbf{P} \right ] + \gamma^3 ( \kappa_{0F} - \kappa_{0E}) \left [ 2 \mathbf{v} \cdot \mathbf{F}/c^3 - ( \mathbf{vv}) : \mathbf{P}/c^3 \right ] \mathbf{v}

where

kappa_{0P	comoving-frame Planck function weighted opacity
kappa_{0E	comoving-frame radiation energy weighted opacity
kappa_{0F	comoving-frame radiation flux weighted opacity

• If the radiation has a blackbody spectrum then kappa_{0E}=\kappa_{0P
• If the radiation is optically thick, then

mathbf{F_0}(\nu_0) \propto -\nabla E_0(\nu_0)/\kappa_0(\nu_0) \propto -[\partial B(\nu_0, T_0) / \partial T_0](\nabla T_0)/\kappa_0(\nu_0

which implies that

kappa_{0F}^{-1}=\kappa_{0R}^{-1}=\frac{\int_0^\infty d \nu_0 \kappa_0(\nu_0)^{-1}[\partial B(\nu_0,T_0)/\partial T_0]}{\int_0^\infty d \nu_0 [\partial B(\nu_0, T_0)/\partial T_0]

• In the optically thin regime, \mathbf{F}_0(\nu_0)| \rightarrow cE_0(\nu_0 so we would have kappa_{0F}=\kappa_{0E however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that kappa_{0F}=\kappa_{0R

Flux limited diffusion

The flux limited diffusion approximation drops the radiation momentum equation in favor of mathbf{F}_0=-\frac{c\lambda}{\kappa_{0R}}\nabla E_

where lambd is the flux-limiter

lambda = \frac{1}{R} \left ( \coth R - \frac{1}{R} \right )

=\frac{|\nabla E_0|}{\kappa_{0R}E_0

which corresponds to a pressure tensor

mathbf{P}_0=\frac{E_0}{2}[(1-R_2)\mathbf{I}+(3R_2-1)\mathbf{n}_0\mathbf{n}_0

_2=\lambda+\lambda^2R^

If we Lorentz boost the comoving terms into the lab frame and keep terms necessary to maintain accuracy we get:

^0=\kappa_{0P} \left ( E-\frac{4 \pi B}{c} \right ) + \left ( \frac{\lambda}{c} \right ) \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \mathbf{v} \cdot \nabla E - \frac{\kappa_{0P}}{c^2} E \left [ \frac{3-R_2}{2}v^2 + \frac{3R_2-1}{2}(\mathbf{v} \cdot \mathbf{n})^2 \right ] + \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \kappa_{0P} \left ( E - \frac{4 \pi B}{c} \right )

mathbf{G} = -\lambda \nabla E + \kappa_{0P} \frac{\mathbf{v}}{c} \left ( E - \frac{4 \pi B}{c} \right ) - \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \lambda \nabla E + 2 \lambda \left ( \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \frac{(\mathbf{v} \cdot \nabla E )\mathbf{v}}{c^2

Numerics of Flux Limited Diffusion

If we plug the expressions for the radiation 4-momentum back into the gas equations and keep terms necessary to maintain accuracy we get:

frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P\color{green}{-\lambda \nabla E

frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = \color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{green}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E} \color{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E

frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} = \color{red}{\kappa_{0P} (4 \pi B-cE)} \color{green}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )} \color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}

For static diffusion, the terms in blue with v²/c can be dropped and the system can be split into the usual hydro update (black), radiative source terms (green) using time centered radiation energy, and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature). If the opacity is independent of temperature and radiation energy density, then the implicit solve only involves the radiation energy density. Otherwise some sort of sub-cycling would be required.

Operator splitting

Krumholz et al. perform Implicit Radiative, Explicit Hydro, Explicit Radiative In AstroBEAR this would look like:

• Initialization
  • Prolongate, d, p, e, E, Edot
• Step 1
  • Overlap d, p, e, E and do physical BC's
  • Do IR which updates e₀, and E₀ using d₁, e₁, E₁
  • Update E_2*mbc using Edot_2*mbc
  • Update e_2*mbc using E_2*mbc, Edot_2*mbc, and e_2*mbc
  • Update Edot₀ using pre IR and post IR E₀
  • Ghost e_2*mbc, E_mbc+1, Edot_mbc+1
  • Do first EH_mbc
  • Do ER_mbc —- Terms with grad E can be done without ghosting since EH did not change E. The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.
  • Store Edot in child arrays to be prolongated
• Step 2
  • Overlap d, p, e, E, and do physical BCs
  • Do IR which updates e₀, and E₀ using d₁, e₁, E₁
  • Update Edot₀ using pre IR and post IR E₀
  • Update E₁ using Edot₁
  • Ghost e_mbc, E₁, Edot₁
  • Do second EH₀
  • Do ER₀ —- Terms with grad E can be done without ghosting since EH did not change E. The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.

Implicit Discretization

For now we will assume that kappa_{0P and kappa_{0R are constant over the implicit update. In this case we can solve the radiation energy equation:

frac{\partial E}{\partial t} = \nabla \cdot \mathbf{F} + \kappa_{0P} (4 \pi B-cE) = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B-cE

frac{\partial e}{\partial t} = - \kappa_{0P} (4 \pi B-cE

where mathbf{F} = \frac{c\lambda}{\kappa_{0R}} \nabla

1D

Which we can discretize for (1D) as

^{n+1}_i-E^{n}_i = \left [ \alpha^n_{i+1/2} \left ( E^{n+1}_{i+1}-E^{n+1}_{i} \right ) - \alpha^n_{i-1/2} \left ( E^{n+1}_{i}-E^{n+1}_{i-1} \right ) \right ] + \epsilon^n_i \left ( \frac{4 \pi}{c} B(T^n_i)-E^{n+1}_i \right )

^{n+1}_i-e^{n}_i = - \epsilon^n_i \left ( \frac{4 \pi}{c} B(T^n_i)-E^{n+1}_i \right )

where

frac{\Delta t}{\Delta x}\mathbf{F}^n_{i+1/2} = \alpha^n_{i+1/2} \left ( E^{n+1}_{i+1} - E^{n+1}_i \right )

and

epsilon^n_i=c\Delta t \kappa^n_{0P,i

represents the number of absorption/emissions during the time step

and the diffusion coefficient is given by

alpha_{i+1/2}=\frac{\Delta t}{\Delta x^2} \frac{c \lambda_{i+1/2}}{\kappa_{i+1/2}} \mbox{ where } \kappa_{i+1/2} = \frac{\kappa_{i}+\kappa_{i+1}}{2} \mbox{ and } \lambda_{i+1/2} = f \left ( R_{i+1/2} \right

where

_{i+1/2} = \frac{\left | E_{i+1}-E_{i} \right | }{2 \kappa_{i+1/2} \left ( E_i+E_{i+1} \right )

This gives matrix coefficients

left ( 1 + \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) E^{n+1}_i - \left ( \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1}=E^n_i+\frac{4\pi \epsilon^n_i}{c}B \left (T^n_i \right )

left ( 1 \right ) e^{n+1}_i - \left ( \epsilon^n_i \right )E^{n+1}_i =e^n_i-\frac{4\pi \epsilon^n_i}{c}B \left (T^n_i \right )

Clearly the second equation is trivial to solve after the first system of equations has been solved. So if we treat the temperature as being constant we can calculate the local heating/cooling due to radiative emission/absorption

^{n+1}_i = e^n_i + \epsilon^n_i \left ( E^{n+1}_i -\frac{4\pi}{c}B \left (T^n_i \right ) \right

Expanding about e₀

Of course the above system is only linear if we ignore the changes in the Temperature due to heating during the implicit step which would feed back into the source function. We can improve this by writing (T) = \frac{c}{4 \pi} a_RT^4 = \frac{c}{4 \pi} \left ( \frac{2}{n f k_B} \right ) ^4 a_R e^4 = \frac{c}{4 \pi} A_R e^4 = A_R e_0^4 \left ( 1 + 4 \delta + 6 \delta^2 + ... \right

where = e_0 \left ( 1 + \delta \right )

Then the system of equations becomes

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + c \kappa_{0P} \left [ A_R e_0^4 \left ( 1 + 4 \delta \right ) -E \right ]

frac{\partial \delta}{\partial t} = - \frac{c \kappa_{0P}}{e_0} \left [ \left (1+4 \delta \right )-E \right ]

or equivalently

frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} \left [ 4 \pi B \left ( T_0 \right ) \left ( 1 + 4 \delta \right ) -cE \right ]

frac{\partial \delta}{\partial t} = - \frac{ \kappa_{0P}}{e_0} \left [ 4 \pi B \left ( T_0 \right ) \left ( 1+4 \delta \right ) -cE \right ]

2D etc…

For 2D the matrix coefficients would be

+\alpha^n_{i+1/2,j}+\alpha^n_{i-1/2,j}+\alpha^n_{i,j+1/2}+\alpha^n_{i,j-1/2}+\epsilon^n_i	^{n+1}_{i,j
\alpha^n_{i-1/2,j	^{n+1}_{i-1,j
\alpha^n_{i+1/2,j	^{n+1}_{i+1,j
\alpha^n_{i,j-1/2	^{n+1}_{i,j-1
\alpha^n_{i,j+1/2	^{n+1}_{i,j+1

and the source term would be unchanged from 1D

^n_i+\frac{4\pi \epsilon^n_i}{c}B \left (T^n_i \right )

Coarse-Fine Boundaries

Since we are doing our implicit solves first, we can use time interpolated solutions for the implicit solve for non-refined ghost zones. To do this we just need Edot. The opacities etc… in the ghost zones can be obtained from the hydro terms.

What about radiative implicit heating in coarse ghost cells which we will need to do hydro update? Could do this explicitly using E and Edot

Physical Boundary Conditions

Open boundaries

We would like the radiation to leave at the free streaming limit. So frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x

Clearly if we set _g = and lambda_g =\kappa_g \Delta we should get the correct flux.

This corresponds to an alpha = c \frac{\Delta t}{\Delta x

So we would just modify alph and zero out the matrix coefficient to the ghost zone

Thermally emitting boundary

Another possible boundary condition would be to have the edge of the grid be adjacent to some thermal emitter. This corresponds to setting the radiation energy density in the ghost zones to the Planck function.

_g = \frac{4 \pi}{c} B(T_g however, we still need an opacity which could be defined from the fluid properties of the ghost region… This would essentially be a thermally emitting boundary with the temperature, density, opacity, etc… derived from the hydro boundary type. If the boundary type was extrapolating, and the density and temperature uniform, one could have a constant thermal energy spectrum… However we would need to either specify the opacity and temperature of the ghost zone - or the various fluid properties needed to reconstruct the opacity and temperature of the ghost zone - or make an extra call to set physicalBC - since this may happen in between a hydro step and a physical boundary update.

Reflecting Boundary

Reflecting boundary should be fairly straightforward. This an be achieved by setting alpha_g = which zeros out any flux - and has the same effect as setting ^{n+1}_g=E^{n+1}_

Constant radiative flux

To have a constant radiative flux we must have alpha_g \left ( E^{n+1}_i-E^{n+1}_g \right ) = F_0 \frac{\Delta t}{\Delta x

Which we can solve for ^{n+1}_g = E^{n+1}_i - \frac{F_0 \Delta t}{\alpha_g \Delta x

but when we plug this into the coefficient matrix the terms with alpha_ cancel and we just get _0 \frac{\Delta t}{\Delta x in the source vector

Summary

Boundary	alph
Open	\frac{\Delta t}{\Delta x	0
User-Defined opacity and Temperature	alpha_	alpha_g \frac{4 \pi}{c} B(T_g
Reflecting
User-Defined Flux		_0 \frac{\Delta t}{\Delta x