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.

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 \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B-cE

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 )

where

epsilon^n_i=\frac{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 )

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 )