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

which if we plug 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₁ using Edot₁
  • 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.

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

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 )

and 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

• At coarse fine boundaries, the simplest approach is to use time interpolated energy fields and opacities in the ghost zones.
• Or we could use time interpolated - or averaged fluxes. Then we don't need fluid variables rho, \dot{\rho}, E, \dot{E to get opacities and boundary fields… We just need the normal fluxes
  • Unfortunately during the implicit solve, there is no way to back out the radiative fluxes directly. So we would have to approximate them - and we need the energy and opacity for the explicit updates anyways…

Hydro before implicit

• First hydro update using extended ghost zones
  • Just need to update ghost opacities internally and in first row of ghost zones
  • Do implicit solve
  • Ghost
• Second hydro update
  • Need to update opacities internally and in first row of ghost zones and then communicate ghost opacities.
  • Need to update radiation field in ghost zones
• Need to update and communicate ghost opacities at internal boundaries

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

So in 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