An Implicit Finite Volume Scheme to Solve the Time-dependent Radiation Transport Equation Based on Discrete Ordinates

We describe the radiation field using frequency-integrated specific intensity I in the lab frame, which is a function of time t, spatial location (x, y, z), and angular direction as defined by the unit vector n . The transport term is greatly simplified in the lab frame, while the source terms describing the interactions between radiation and gas are convenient in the comoving frame of the fluid. Therefore, we also adopt this mixed frame approach (Lowrie et al. 1999; Jiang et al. 2012, 2014). The RT equation we solve is (Mihalas & Klein 1982; Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&\displaystyle \frac{\partial I}{\partial t}+c{\boldsymbol{n}}\cdot {\boldsymbol{\nabla }}I=c\left(\eta -\chi I\right),\end{eqnarray} \tag{ 1 }$$

where c is the speed of light. The emissivity η and opacity χ are formally defined in the lab frame in the above equation. But the whole source term is calculated in the comoving frame via the approach below. Traditionally, in order to get the expression for the source term, expansion to the first order of ${ \mathcal O }\left(v/c\right)$ is usually used (Mihalas & Klein 1982), with some additional second-order terms (Lowrie et al. 1999; Krumholz et al. 2007; Jiang et al. 2014). Here we apply the Lorentz transformation to the specific intensities to completely avoid the need for expansion in terms of v/c. Although this may not sound necessary for nonrelativistic systems, this approach does not cause additional computational cost, and it results in a much cleaner calculation of the source terms.

The frequency-integrated specific intensities in the lab frame I( n ) are related to the comoving frame values ${I}_{0}({{\boldsymbol{n}}}^{{\prime} })$ via the Lorentz transformation as (Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&{I}_{0}({{\boldsymbol{n}}}^{{\prime} })={\gamma }^{4}{\left(1-{\boldsymbol{n}}\cdot {\boldsymbol{v}}/c\right)}^{4}I({\boldsymbol{n}})\equiv {{\rm{\Gamma }}}^{4}({\boldsymbol{n}},{\boldsymbol{v}})I({\boldsymbol{n}}),\end{eqnarray} \tag{ 2 }$$

where v is the flow velocity, $\gamma \equiv 1/\sqrt{1-{v}^{2}/{c}^{2}}$ is the corresponding Lorentz factor, and ${{\boldsymbol{n}}}^{{\prime} }$ is the angle in the comoving frame given by

$$\begin{eqnarray}&&{{\boldsymbol{n}}}^{{\prime} }=\displaystyle \frac{1}{\gamma \left(1-{\boldsymbol{n}}\cdot {\boldsymbol{v}}/c\right)}\left[{\boldsymbol{n}}-\gamma \displaystyle \frac{{\boldsymbol{v}}}{c}\left(1-\displaystyle \frac{\gamma }{\gamma +1}\displaystyle \frac{{\boldsymbol{n}}\cdot {\boldsymbol{v}}}{c}\right)\right].\end{eqnarray} \tag{ 3 }$$

We have defined ${\rm{\Gamma }}({\boldsymbol{n}},{\boldsymbol{v}})\equiv \gamma \left(1-{\boldsymbol{n}}\cdot {\boldsymbol{v}}/c\right)$ for convenience. We multiply the left-hand sides of the source term with Γ⁴ to get

$$\begin{eqnarray}&&\displaystyle \frac{\partial {I}_{0}}{\partial t}=c{\rm{\Gamma }}\left[{\eta }_{0}-{\chi }_{0}{I}_{0}\right].\end{eqnarray} \tag{ 4 }$$

Here the frequency-integrated emissivity (η₀) and opacity (χ₀) in the comoving frame are related to η, χ as

$$\begin{eqnarray}&&{\eta }_{0}={{\rm{\Gamma }}}^{3}\eta ,\,{\chi }_{0}={{\rm{\Gamma }}}^{-1}\chi .\end{eqnarray} \tag{ 5 }$$

Under the assumption of local thermal equilibrium, we can write the frequency-integrated source term in the comoving frame as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {I}_{0}}{\partial t} & = & c{\rm{\Gamma }}\left[\rho {\kappa }_{s}\left({J}_{0}-{I}_{0}\right)\right.\\ & & +\ \left.\rho {\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{T}^{4}}{4\pi }-{I}_{0}\right)+\rho {\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{T}^{4}}{4\pi }-{J}_{0}\right)\right]\\ & = & \ c{\rm{\Gamma }}\left[\rho ({\kappa }_{s}+{\kappa }_{a})\left({J}_{0}-{I}_{0}\right)\right.\\ & & +\ \left.\rho \left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left(\displaystyle \frac{{a}_{r}{T}^{4}}{4\pi }-{J}_{0}\right)\right]\end{array}\end{eqnarray} \tag{ 6 }$$

and the corresponding equation for gas internal energy E_g as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{g}}{\partial t}=-c\rho \left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left({a}_{r}{T}^{4}-4\pi {J}_{0}\right).\end{eqnarray} \tag{ 7 }$$

Here a_r is the radiation constant, and ρ and T are the gas density and temperature, which should be comoving frame quantities, in principle. When the RT algorithm is coupled to a Newtonian MHD solver, we do not distinguish gas quantities as well as time in the lab and comoving frames. The angular averaged mean intensity in the comoving frame J₀ is defined as

$$\begin{eqnarray}&&{J}_{0}\equiv \displaystyle \frac{1}{4\pi }\int {I}_{0}d{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 8 }$$

where Ω₀ is the angular element in the comoving frame. The specific scattering opacity κ_s (with units of cm² g⁻¹) is assumed to be isotropic and does not change the mean photon energy such as electron scattering. The Rosseland mean absorption opacity is κ_a , and we take κ_{δ P} to be the difference between the Planck mean and Rosseland mean opacity. This is designed to ensure that when we integrate the zeroth and first moment of the comoving frame specific intensity over the angular space, we get the energy and momentum source terms $\rho \left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left({a}_{r}{T}^{4}-4\pi {J}_{0}\right)$ and $-\rho \left({\kappa }_{s}+{\kappa }_{a}\right)\left(4\pi {{\boldsymbol{H}}}_{0}\right)$ with ${{\boldsymbol{H}}}_{0}\equiv \int {I}_{0}{{\boldsymbol{n}}}^{{\prime} }d{{\rm{\Omega }}}_{0}/\left(4\pi \right)$. This is consistent with the Rosseland and Planck mean opacities normally used for radiation moment equations. We multiply both sides of Equation (6) with Γ⁻⁴ to get the source term for lab frame specific intensities.

Photons are coupled to the gas via momentum and energy exchange. Since the total energy and momentum are conserved in the lab frame, we integrate the source terms for lab frame specific intensities over the angles to get the full radiation MHD equations as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\\ \displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{v}}{\boldsymbol{v}}-{\boldsymbol{B}}{\boldsymbol{B}}+{{\mathsf{P}}}^{* }\right)=-{{\boldsymbol{S}}}_{{\boldsymbol{r}}}({\boldsymbol{P}}),\\ \displaystyle \frac{\partial E}{\partial t}+{\boldsymbol{\nabla }}\cdot \left[\left(E+{P}^{* }\right){\boldsymbol{v}}-{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{v}})\right]=-{S}_{r}(E),\\ \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}-{\boldsymbol{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})=0,\\ \displaystyle \frac{\partial I}{\partial t}+c{\boldsymbol{n}}\cdot {\boldsymbol{\nabla }}I={{cS}}_{I},\\ {S}_{I}\equiv {{\rm{\Gamma }}}^{-3}\left[\rho \left({\kappa }_{s}+{\kappa }_{a}\right)\left({J}_{0}-{I}_{0}\right)\right.\\ \ +\ \left.\rho ({\kappa }_{a}+{\kappa }_{\delta P})\left(\displaystyle \frac{{a}_{r}{T}^{4}}{4\pi }-{J}_{0}\right)\right],\\ {S}_{r}(E)\equiv 4\pi c\displaystyle \int {S}_{I}d{\rm{\Omega }},\\ {{\boldsymbol{S}}}_{{\boldsymbol{r}}}({\boldsymbol{P}})\equiv 4\pi \displaystyle \int {\boldsymbol{n}}{S}_{I}d{\rm{\Omega }}.\end{array}\end{eqnarray} \tag{ 9 }$$

Here B is the magnetic field, P ^* ≡ (P + B²/2)I (with I the unit tensor), and the magnetic permeability μ is taken to be 1 in this unit system. The total gas energy density is

$$\begin{eqnarray}&&E={E}_{g}+\displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{{B}^{2}}{2}.\end{eqnarray} \tag{ 10 }$$

We assume ideal gas so that the gas internal energy is related to gas pressure via the adiabatic index γ_g as E_g = P/(γ_g − 1) for γ_g ≠ 1. The gas temperature is calculated with $T=P/\left({R}_{\mathrm{ideal}}\rho \right)$, where R_ideal is the ideal gas constant.

The radiation field shares the same spatial grid as used for the gas quantities in Athena++. We label the three coordinates as x, y, z, which can be any orthogonal curvilinear systems such as Cartesian, cylindrical, or spherical polar coordinates. The total number of grid points along the three directions are N_x , N_y , and N_z , respectively. Notice that our finite volume scheme does not need the grid spacing to be uniform. What we need is the volume of each cell V_i,j,k for i ∈ [0, N_x − 1], j ∈ [0, N_y − 1], and k ∈ [0, N_z − 1], as well as surface areas of each cell facing the three directions _Ax , A_y , and A_z . All of the specific intensities are defined at the volume center of each cell.

We adopt the same angular discretization scheme as used by Davis et al. (2012) and Jiang et al. (2014) to be our default choice. ¹ The angular unit vectors n are defined with respect to a global Cartesian coordinate, and they do not change with spatial locations. The three components of each angle ${\boldsymbol{n}}=\left({\mu }_{x},{\mu }_{y},{\mu }_{z}\right)$, as well as the quadrature weight w_n , are determined based on the same algorithm developed by Bruls et al. (1999) following the original quadrature principle described in Carlson (1963). This set of angles ensures that the following relations are satisfied numerically:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{n=0}^{N-1}{w}_{n} & = & 1,\\ \displaystyle \sum _{n=0}^{N-1}{\mu }_{x}{w}_{n} & = & \displaystyle \sum _{n=0}^{N-1}{\mu }_{y}{w}_{n}=\displaystyle \sum _{n=0}^{N-1}{\mu }_{z}{w}_{n}=0,\\ \displaystyle \sum _{n=0}^{N-1}{\mu }_{x}^{2}{w}_{n} & = & \displaystyle \sum _{n=0}^{N-1}{\mu }_{y}^{2}{w}_{n}=\displaystyle \sum _{n=0}^{N-1}{\mu }_{z}^{2}{w}_{n}=\displaystyle \frac{1}{3},\end{array}\end{eqnarray} \tag{ 11 }$$

where N is the total number of angles.

Notice that these angles are always defined in the lab frame. The angles in the comoving frame are determined based on Equation (3), while the quadrature weight in the comoving frame is ${w}_{n}^{{\prime} }={{\rm{\Gamma }}}^{-2}{w}_{n}$. For an infinite number of angles, it can be shown that ∫Γ⁻² dΩ₀ = 1 for any flow velocity. However, this relation normally is not satisfied numerically to machine precision due to truncation errors with a finite number of angles. We correct this error by normalizing ${w}_{n}^{{\prime} }$ as

$$\begin{eqnarray}&&{w}_{n}^{{\prime} }\equiv \displaystyle \frac{{{\rm{\Gamma }}}^{-2}{w}_{n}}{{\displaystyle \sum }_{n=0}^{N-1}{{\rm{\Gamma }}}^{-2}{w}_{n}}.\end{eqnarray} \tag{ 12 }$$

Therefore, ${\sum }_{0}^{N-1}{w}_{n}^{{\prime} }=1$ is always satisfied. We find it is necessary to correct this angular truncation error to avoid an accumulative energy error when the flow velocity is not zero. The commonly used radiation energy density, radiation flux, and radiation pressure in the lab and comoving frames are just derived quantities from specific intensities as

$$\begin{eqnarray}\begin{array}{rcl}{E}_{r} & = & 4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{n}{w}_{n},{E}_{r,0}=4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{0,n}{w}_{n}^{{\prime} },\\ {{\boldsymbol{F}}}_{r} & = & 4\pi c\displaystyle \sum _{n=0}^{N-1}{I}_{n}{\boldsymbol{n}}{w}_{n},{{\boldsymbol{F}}}_{r,0}=4\pi c\displaystyle \sum _{n=0}^{N-1}{I}_{0,n}{{\boldsymbol{n}}}^{{\prime} }{w}_{n}^{{\prime} },\\ {{\mathsf{P}}}_{r} & = & 4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{n}{\boldsymbol{n}}{\boldsymbol{n}}{w}_{n},{{\mathsf{P}}}_{r,0}=4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{0,n}{{\boldsymbol{n}}}^{{\prime} }{{\boldsymbol{n}}}^{{\prime} }{w}_{n}^{{\prime} }.\end{array}\end{eqnarray} \tag{ 13 }$$

The radiation energy source term can be very stiff in the optically thick regions, while the radiation momentum source term is typically not (Sekora & Stone 2010; Jiang et al. 2012). Therefore, when we solve the RT equation, we will assume that the flow velocity and gas density do not change. But we evolve the gas temperature with the radiation variables together to ensure that correct thermal equilibrium state is achieved. As higher-order spatial reconstruction without a flux limiter cannot guarantee total variation diminishing, which is necessary to avoid oscillatory behavior around local maxima or minima, we will only use first-order spatial reconstruction for our algorithm. The commonly used flux limiters will make the transport term nonlinear, which is very hard to solve implicitly. The algorithm will also be first-order accurate in time for RT to avoid oscillatory behavior during the temporal evolution of gas temperature and radiation energy density when the source term is very stiff.

For the given initial conditions ${I}_{n}^{m},{\rho }^{m},{T}^{m},{v}^{m}$ at time step m, the specific intensities ${I}_{n}^{m+1}$ with n ∈ [0, N − 1] at time step m + 1 are calculated based on the following N + 1 equations:

$$\begin{eqnarray}&&\begin{array}{l}{I}_{n}^{m+1}-{I}_{n}^{m}+{\rm{\Delta }}{tc}{\boldsymbol{n}}\cdot {\boldsymbol{\nabla }}{I}_{n}^{m+1}={\rm{\Delta }}{tc}{{\rm{\Gamma }}}_{n}^{-3}\\ \ \times \ \left[{\rho }^{m}{\kappa }_{s}\left({J}_{0}^{m+1}-{I}_{0,n}^{m+1}\right)\right.\\ \ \ +\ {\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}^{m+1}\right)}^{4}}{4\pi }-{I}_{0,n}^{m+1}\right)\\ \ \ +\ \left.{\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}^{m+1}\right)}^{4}}{4\pi }-{J}_{0}^{m+1}\right)\right],\\ \ \ \times \ \displaystyle \frac{{\rho }^{m}{R}_{\mathrm{ideal}}}{{\gamma }_{g}-1}\left({T}^{m+1}-{T}^{m}\right)\\ \ =\ -{\rm{\Delta }}{tc}{\rho }^{m}\left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left[{a}_{r}{\left({T}^{m+1}\right)}^{4}-4\pi {J}_{0}^{m+1}\right],\end{array}\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{I}_{0,n}^{m+1}={{\rm{\Gamma }}}_{n}^{4}{I}_{n}^{m+1},{J}_{0}^{m+1}=\displaystyle \sum _{n=0}^{N-1}{w}_{n}^{{\prime} }{I}_{0,n}^{m+1},\end{eqnarray} \tag{ 15 }$$

and Γ_n is calculated using v^m for each angle. The opacities (κ_s , κ_a , and κ_{δ P} ) can be functions of local gas properties, in principle. But we will assume that they keep the values at time step m in the above equations. This is a set of fully coupled equations for gas temperature and specific intensities at all angles (via the terms with ${J}_{0}^{m+1}$) and spatial locations (via the transport term ${\boldsymbol{n}}\cdot {\boldsymbol{\nabla }}{I}_{n}^{m+1}$). These equations are only linear with respect to ${I}_{n}^{m+1}$ but not for T^m+1 due to the ${({T}^{m+1})}^{4}$ term. All of the transport and source terms are calculated using the variables at time step m + 1 to make the scheme unconditionally stable for any time step Δt. We will use an iterative approach to solve these equations as described below.

3.2.1. The Transport Term

Following Jiang et al. (2014), we can rewrite the transport term in the above equations as $c{\boldsymbol{\nabla }}\cdot \left({\boldsymbol{n}}{I}_{n}^{m+1}\right)$ with our default angular discretization because n does not change with spatial location. Then this term becomes the conservative format and can be calculated with the standard finite volume approach. For the interface at i − 1/2, j, k, we simply take the left and right states for ${I}_{n}^{m+1}$ to be ${I}_{n}^{m+1}(i-1,j,k)$ and ${I}_{n}^{m+1}(i,j,k)$, respectively. The same approach is used for other directions.

Special care is needed to determine the surface fluxes for each specific intensity. If we simply take the upwind fluxes as in Stone & Mihalas (1992), although the total energy is conserved, numerical diffusion can overcome the physical diffusion in the very optically thick regime. This can be easily demonstrated with only two angles in one dimension. If we assume the two angles have μ_x,1 > 0 and μ_x,2 < 0, the upwind fluxes at i − 1/2 will be c μ_x,1 I₁(i − 1) and c μ_x,2 I₂(i). Similarly, the fluxes at i + 1/2 for the two angles are c μ_x,1 I₁(i) and c μ_x,2 I₂(i + 1). The change of the mean radiation energy density w₁ I₁ + w₂ I₂ at cell i is then determined by cw₂ μ_x,2 I₂(i + 1) + cw₁ μ_x,1 I₁(i) − cw₂ μ_x,2 I₂(i) − cw₁ μ_x,1 I₁(i − 1). However, the radiation moment equation requires that it should be proportional to F_r (i + 1) − F_r (i − 1) = cw₂ μ_x,2 I₂(i + 1) + cw₁ μ_x,1 I₁(i + 1) − cw₂ μ_x,2 I₂(i − 1) − cw₁ μ_x,1 I₁(i − 1). The difference is essentially the truncation error, which is negligible in the optically thin regime. However, in the optically thick regime, when the radiation field is very close to being isotropic, ∣I₂(i + 1) − I₁(i + 1)∣ can be much smaller than ∣I₁(i + 1) − I₁(i)∣. This will cause the numerical flux due to the truncation errors to be much larger than the physical flux due to photon diffusion, which results in completely wrong numerical solutions. We have designed a modified HLLE flux to overcome the issue.

The dynamical diffusion regime also needs to be treated carefully, as discussed in Jiang et al. (2014). The numerical truncation error due to finite flow velocity can also overcome the physical diffusion when the optical depth per cell is very large, particularly with first-order spatial reconstruction. In principle, these numerical errors can be minimized with second- or third-order spatial reconstructions. However, it is not an option for our implicit scheme. Instead, we find that a specially designed HLLE-type flux can work very well in all regimes.

We rewrite the transport term as

$$\begin{eqnarray}&&\begin{array}{l}c{\boldsymbol{\nabla }}\cdot \left({\boldsymbol{n}}{I}_{n}^{m+1}\right)={\boldsymbol{\nabla }}\cdot \left[\left(c{\boldsymbol{n}}-f{{\boldsymbol{v}}}^{m}\right){I}_{n}^{m+1}\right]\\ \qquad \ +\ {\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m+1}\right).\end{array}\end{eqnarray} \tag{ 16 }$$

Here f is a function of optical depth per cell ${\tau }_{c}\,\equiv {\rho }^{m}\left({\kappa }_{a}+{\kappa }_{s}\right){\rm{\Delta }}L$, with ΔL to be a constant factor times the cell size. The function is chosen as

$$\begin{eqnarray}&&f=1-\exp (-{\tau }_{c}^{2}).\end{eqnarray} \tag{ 17 }$$

This is designed to separate the advection term when the photon mean free path is not resolved (f ≈ 1). This is not necessary when the optical depth per cell is much smaller than 1 (f ≈ 0). Since our time step will satisfy the CFL condition for flow velocity, we can calculate the advective term ${\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}^{m+1}\right)$ explicitly, which means

$$\begin{eqnarray}&&\begin{array}{l}c{\boldsymbol{\nabla }}\cdot \left({\boldsymbol{n}}{I}_{n}^{m+1}\right)\approx {\boldsymbol{\nabla }}\cdot \left[\left(c{\boldsymbol{n}}-f{{\boldsymbol{v}}}^{m}\right){I}_{n}^{m+1}\right]\\ \qquad \ +\ {\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right).\end{array}\end{eqnarray} \tag{ 18 }$$

We use second-order spatial reconstruction for ${I}_{n}^{m}$ to get the values at the cell faces and then upwind flux based on f v ^m to calculate the term ${\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right)$.

For the interface at (i − 1/2, j, k), we take the left and right states for each specific intensity to be ${I}_{n,L}={I}_{n}^{m+1}(i-1,j,k)$ and ${I}_{n,R}={I}_{n}^{m+1}(i,j,k)$. Then we calculate the flux F_n (i − 1/2) at the interface for each angle based on a modified HLLE formula (Appendix). The key modification is to make sure the flux is the upwind value in the optically thin case and gets close to the centered value in the optically thick case based on the parameter τ_c . Similarly, we calculate the flux F_n (i + 1/2) for the interface (i + 1/2, j, k), as well as fluxes along other directions F_n (j − 1/2), F_n (j + 1/2), F_n (k − 1/2), and F_n (k + 1/2). Detailed expressions for the fluxes are given in the Appendix. Then the flux divergence term is calculated as

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{\nabla }}\cdot \left[\left(c{\boldsymbol{n}}-f{{\boldsymbol{v}}}^{m}\right){I}_{n}^{m+1}\right]=\displaystyle \frac{1}{{V}_{i,j,k}}\\ \ \times \ \left[\right.{A}_{x}(i+1/2){F}_{n}(i+1/2)-{A}_{x}(i-1/2){F}_{n}(i-1/2)\\ \ +\ {A}_{y}(j+1/2){F}_{n}(j+1/2)-{A}_{y}(j-1/2){F}_{n}(j-1/2)\\ \ +\ \left.{A}_{z}(k+1/2){F}_{n}(k+1/2)-{A}_{z}(k-1/2){F}_{n}(k-1/2)\right].\end{array}\end{eqnarray} \tag{ 19 }$$

Reorganizing all terms results in the following equation that we will apply our iterative solver:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{g}_{1}{I}_{n}^{m+1}+{g}_{2}{I}_{n}^{m+1}(i-1)+{g}_{3}{I}_{n}^{m+1}(i+1)\\ \,+\,{g}_{4}{I}_{n}^{m+1}(j-1)+{g}_{5}{I}_{n}^{m+1}(j+1)\\ \,+\,{g}_{6}{I}_{n}^{m+1}(k-1)+{g}_{7}{I}_{n}^{m+1}(k+1)\\ \,=\,{I}_{n}^{m}-{\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right)\\ \,+\,{\rm{\Delta }}tc{{\rm{\Gamma }}}_{n}^{-3}\left[\right.{\rho }^{m}{\kappa }_{s}\left({J}_{0}^{m+1}-{I}_{0,n}^{m+1}\right)\\ \,+\,{\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}^{m+1}\right)}^{4}}{4\pi }-{I}_{0,n}^{m+1}\right)\\ \,+\,\left.{\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}^{m+1}\right)}^{4}}{4\pi }-{J}_{0}^{m+1}\right)\right],\\ \,\displaystyle \frac{{\rho }^{m}{R}_{{\rm{i}}{\rm{d}}{\rm{e}}{\rm{a}}{\rm{l}}}}{{\gamma }_{g}-1}\left({T}^{m+1}-{T}^{m}\right)\\ \,=\,-{\rm{\Delta }}tc{\rho }^{m}\left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left[{a}_{r}{\left({T}^{m+1}\right)}^{4}-4\pi {J}_{0}^{m+1}\right].\end{array}\end{array}\end{eqnarray} \tag{ 20 }$$

Expressions for the coefficients g₁–g₇ are given in the Appendix. Notice that their values depend on the local optical depth τ_c , and they only need to be calculated once for each time step.

3.2.2. The Iterative Solver

We need to solve the N + 1 Equation (20) simultaneously to get the solutions for ${I}_{n}^{m+1}$ (n ∈ [0, N − 1]) and T^m+1. These equations are only linear for ${I}_{n}^{m+1}$. These terms can be rewritten as a sparse matrix with a size of N × N_x × N_y × N_z . If the specific intensities are ordered in a 1D array such that ${I}_{n}^{m+1}(i,j,k)$ is located at kN_y N_x N + jN_x N + iN + n, the matrix has the special structure that the diagonal direction are blocks with size N × N, which couple all of the angles at the same spatial locations. All other elements are only nonzero at locations corresponding to (n, i ± 1, j, k), (n, i, j ± 1, k), and (n, i, j, k ± 1). The terms related to T^m+1 are nonlinear, but T^m+1 at different spatial locations are completely decoupled.

This motivates us to solve these equations iteratively, similar to the Jacobi approach used for linear systems. For each cell, we take all of the specific intensities at different spatial locations to be the values from the last step. Then we go through each cell and solve the equations for all of the specific intensities at once. We update the off-diagonal terms with the updated specific intensities and continue this process until the changes of the specific intensities between neighboring steps are below a certain threshold.

Specifically, we solve the following equations iteratively:

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{g}_{1}{I}_{n,l}^{m+1}+{g}_{2}{I}_{n,l-1}^{m+1}(i-1)+{g}_{3}{I}_{n,l-1}^{m+1}(i+1)\\ \,+\,{g}_{4}{I}_{n,l-1}^{m+1}(j-1)+{g}_{5}{I}_{n,l-1}^{m+1}(j+1)\\ \,+\,{g}_{6}{I}_{n,l-1}^{m+1}(k-1)+{g}_{7}{I}_{n,l-1}^{m+1}(k+1)\\ ={I}_{n}^{m}-{\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right)\\ \,+\,{\rm{\Delta }}tc{{\rm{\Gamma }}}_{n}^{-3}\left[\right.{\rho }^{m}{\kappa }_{s}\left({J}_{0,l}^{m+1}-{I}_{0,n,l}^{m+1}\right)\\ \,+\,{\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{I}_{0,n,l}^{m+1}\right)\\ \,+\,{\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{J}_{0,l}^{m+1}\right)\left.\right],\\ \,\displaystyle \frac{{\rho }^{m}{R}_{{\rm{i}}{\rm{d}}{\rm{e}}{\rm{a}}{\rm{l}}}}{{\gamma }_{g}-1}\left({T}_{l}^{m+1}-{T}^{m}\right)\\ \,=\,-{\rm{\Delta }}tc{\rho }^{m}\left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left[{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}-4\pi {J}_{0,l}^{m+1}\right],\end{array}\end{array}\end{eqnarray} \tag{ 21 }$$

where l represents each step during the iterative process. At the beginning, we need an initial guess for specific intensities at the whole grid for all of the angles. This is usually taken to be the solution from the last time step or the initial condition. Then we go through all of the grid points to solve the above equation for all of the angles at each cell together in the following way.

For each step during the iteration, since all of the variables ${I}_{n,l-1}^{m+1}(i-1),{I}_{n,l-1}^{m+1}(i+1)$, ${I}_{n,l-1}^{m+1}(j-1),{I}_{n,l-1}^{m+1}(j+1))$, and ${I}_{n,l-1}^{m+1}(k-1),{I}_{n,l-1}^{m+1}(k+1)$ are already known, we can reorganize all of the equations we need to solve at each cell as

$$\begin{eqnarray}\begin{array}{rcl}{g}_{1}{I}_{n,l}^{m+1} & = & {I}_{c}+{\rm{\Delta }}{tc}{{\rm{\Gamma }}}_{n}^{-3}\left[\right.{\rho }^{m}{\kappa }_{s}\left({J}_{0,l}^{m+1}-{I}_{0,n,l}^{m+1}\right)\\ & & \ +\ {\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{I}_{0,n,l}^{m+1}\right)\\ & & \ +\ {\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{J}_{0,l}^{m+1}\right)\left.\right],\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rho }^{m}{R}_{\mathrm{ideal}}}{{\gamma }_{g}-1}\left({T}_{l}^{m+1}-{T}^{m}\right)\\ \ =\ -{\rm{\Delta }}{tc}{\rho }^{m}\left({\kappa }_{a}+{\kappa }_{\delta P}\right)\left[{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}-4\pi {J}_{0,l}^{m+1}\right],\end{array}\end{eqnarray} \tag{ 23 }$$

where I_c includes all other terms that are already known before step l. We multiply the left and right sides of the equation for each ${I}_{n,l}^{m+1}$ with ${{\rm{\Gamma }}}_{n}^{4}$ to get the equation for comoving frame specific intensities as

$$\begin{eqnarray}\begin{array}{rcl}{g}_{1}{I}_{0,n,l}^{m+1} & = & {{\rm{\Gamma }}}_{n}^{4}{I}_{c}+{\rm{\Delta }}{tc}{{\rm{\Gamma }}}_{n}\left[\right.{\rho }^{m}{\kappa }_{s}\left({J}_{0,l}^{m+1}-{I}_{0,n,l}^{m+1}\right)\\ & & \ +\ {\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{I}_{0,n,l}^{m+1}\right)\\ & & \ +\ \left.{\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{J}_{0,l}^{m+1}\right)\right].\end{array}\end{eqnarray} \tag{ 24 }$$

We sum up the equations with the appropriate weight ${w}_{n}^{{\prime} }$ for each angle and get

$$\begin{eqnarray}\begin{array}{rcl}{g}_{1}{J}_{0,l}^{m+1} & = & \displaystyle \sum _{n}{w}_{n}^{{\prime} }{{\rm{\Gamma }}}_{n}^{4}{I}_{c}+{\rm{\Delta }}{tc}{{\rm{\Gamma }}}_{n}\left[\right.\\ & & \ +\ {\rho }^{m}{\kappa }_{a}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{J}_{0,l}^{m+1}\right)\\ & & \ +\ \left.{\rho }^{m}{\kappa }_{\delta P}\left(\displaystyle \frac{{a}_{r}{\left({T}_{l}^{m+1}\right)}^{4}}{4\pi }-{J}_{0,l}^{m+1}\right)\right].\end{array}\end{eqnarray} \tag{ 25 }$$

The scattering term drops out because it does not change the mean photon energy. Together with Equation (23), we can get the solution for ${J}_{0,l}^{m+1}$ and ${T}_{l}^{m+1}$ by performing a root finding for a fourth-order polynomial. Then the specific intensity for each angle ${I}_{0,n,l}^{m+1}$ can be calculated from Equation (24) independently, which can be converted to ${I}_{n,l}^{m+1}$ directly. In this way, the nonlinearity of the source term, as well as the coupling between different angles, all comes down to a fourth-order polynomial solver, which is trivial and independent of the number of angles. Therefore, the whole cost for each step during the iteration is linearly proportional to the number of angles. We continue the process for each grid and stop the iteration when

$$\begin{eqnarray}&&{\rm{\Delta }}I\equiv \displaystyle \sum _{n,i,j,k}| {I}_{n,l}^{m+1}-{I}_{n,l-1}^{m+1}| /\displaystyle \sum _{n,i,j,k}| {I}_{n,l}^{m+1}| \end{eqnarray} \tag{ 26 }$$

is smaller than a threshold, typically 10⁻⁵.

3.2.3. The Radiation Source Terms

Since total energy and momentum are conserved in the lab frame, we use the conservation laws to determine the source terms we need to add to the gas. The radiation energy and momentum at the beginning of the time step are

$$\begin{eqnarray}&&{E}_{r}^{m}=4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{n}^{m}{w}_{n},{{\boldsymbol{F}}}_{r}^{m}=4\pi c\displaystyle \sum _{n=0}^{N-1}{\boldsymbol{n}}{I}_{n}^{m}{w}_{n}.\end{eqnarray} \tag{ 27 }$$

At the end of the time step, they can be calculated as

$$\begin{eqnarray}\begin{array}{rcl}{E}_{r}^{m+1} & = & 4\pi \displaystyle \sum _{n=0}^{N-1}{I}_{n}^{m+1}{w}_{n},\\ {{\boldsymbol{F}}}_{r}^{m+1} & = & 4\pi c\displaystyle \sum _{n=0}^{N-1}{\boldsymbol{n}}{I}_{n}^{m+1}{w}_{n}.\end{array}\end{eqnarray} \tag{ 28 }$$

Part of the change is due to the transport term, which can also be calculated as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{E}_{r} & = & 4\pi \displaystyle \sum _{n=0}^{N-1}\left\{\right.{\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left[\left(c{\boldsymbol{n}}-f{{\boldsymbol{v}}}^{m}\right){I}_{n}^{m+1}\right]\\ & & +\ {\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right)\left.\right\}{w}_{n},\\ {\rm{\Delta }}{{\boldsymbol{F}}}_{r} & = & 4\pi c\displaystyle \sum _{n=0}^{N-1}\left\{\right.{\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left[\left(c{\boldsymbol{n}}-f{{\boldsymbol{v}}}^{m}\right){I}_{n}^{m+1}\right]\\ & & +\ {\rm{\Delta }}t{\boldsymbol{\nabla }}\cdot \left(f{{\boldsymbol{v}}}^{m}{I}_{n}^{m}\right)\left.\right\}{\boldsymbol{n}}{w}_{n}.\end{array}\end{eqnarray} \tag{ 29 }$$

The source terms are then given by

$$\begin{eqnarray}\begin{array}{rcl}{S}_{r}(E) & = & {E}_{r}^{m+1}-{E}_{r}^{m}-{\rm{\Delta }}{E}_{r},\\ {{\boldsymbol{S}}}_{{\boldsymbol{r}}}({\boldsymbol{P}}) & = & \displaystyle \frac{1}{{c}^{2}}\left({{\boldsymbol{F}}}_{r}^{m+1}-{{\boldsymbol{F}}}_{r}^{m}-{\rm{\Delta }}{{\boldsymbol{F}}}_{r}\right).\end{array}\end{eqnarray} \tag{ 30 }$$

3.2.4. Additional Terms with Spherical Polar Angular System

The choice of the angular system is not unique. If the angles vary with spatial locations, additional terms will show up when we convert n · ∇ I to ${\boldsymbol{\nabla }}\cdot \left({\boldsymbol{n}}I\right)$. The RT equation with general angular systems is described in Davis & Gammie (2020). One special case that is useful for nonrelativistic systems is the angular system in spherical polar coordinates, which allows us to use this scheme in 1D or 2D spherical polar coordinates. Notice that our default angular system as described in Section 3.1 can also be used, in principle. The alternative angular system we describe here can be more convenient and efficient for certain applications.

In spherical polar coordinates, we define a set of angles with respect to the local coordinates (r, θ, ϕ) in each cell. We divide the longitudinal angles ψ uniformly from zero to 2π. For the polar angle ζ, it is uniformly divided in $\cos \zeta$ for ζ ∈ [0, π]. These angles point to different directions at different spatial locations, which introduce additional terms in the transport term as (Davis & Gammie 2020)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial I}{\partial t}+c{\boldsymbol{n}}\cdot {\boldsymbol{\nabla }}I=\displaystyle \frac{\partial I}{\partial t}+c{\boldsymbol{\nabla }}\cdot \left({\boldsymbol{n}}I\right)\\ \ -\ \displaystyle \frac{1}{r\sin \zeta }\displaystyle \frac{\partial \left({\sin }^{2}\zeta I\right)}{\partial \zeta }-\displaystyle \frac{\sin \zeta \cos \theta }{r\sin \theta }\displaystyle \frac{\partial \left(\sin \psi I\right)}{\partial \psi }.\end{array}\end{eqnarray} \tag{ 31 }$$

All other terms are unchanged. The numerical scheme described in the previous sections remains the same, except that we need to calculate the change of I due to the new terms.

The new terms represent the transport of specific intensities in the angular space, which we can calculate implicitly as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{I}_{n}^{m+1}-{I}_{n}^{m}}{{\rm{\Delta }}t}-\displaystyle \frac{1}{r}\displaystyle \frac{{\sin }^{2}{\zeta }_{n+1/2}{I}_{n+1/2}^{m+1}-{\sin }^{2}{\zeta }_{n-1/2}{I}_{n-1/2}^{m+1}}{\cos {\zeta }_{n-1/2}-\cos {\zeta }_{n+1/2}}\\ \ -\ \displaystyle \frac{\sin {\zeta }_{n}\cos \theta }{r\sin \theta }\displaystyle \frac{\sin {\psi }_{n+1/2}{I}_{n+1/2}^{m+1}-\sin {\psi }_{n-1/2}{I}_{n-1/2}^{m+1}}{{\psi }_{n+1/2}-{\psi }_{n-1/2}}=0.\end{array}\end{eqnarray} \tag{ 32 }$$

Here the variables with subscripts n − 1/2 and n + 1/2 represent the left- and right-hand sides in the angular space. We use first-order upwind values for ${I}_{n-1/2}^{m+1}$ and ${I}_{n+1/2}^{m+1}$. During each iteration, we modify the coefficient g₁ in Equation (21) to include additional contributions for ${I}_{n,l}^{m+1}$ from this term. For all other specific intensities with different angles in the above equation, we use values from the last iteration, which are added to the right-hand side of Equation (21). We find it necessary to include all of the terms during the iterative process in order to maintain a balance between all of the transport and source terms.

The implicit solver is coupled to the two-step van Leer MHD solver in Athena++ (Stone et al. 2020). We will not repeat the details of the Godunov scheme with constrained transport for MHD here. Instead, we just outline the steps of the full algorithm.

• (a)  
  Determine the time step Δt based on the flow speed, sound speed, Alfvén speed, and cell sizes with CFL number 0.4.
• (b)  
  Update the gas quantities and magnetic field for a half time step Δt/2 with the MHD solver in Athena++.
• (c)  
  Evolve the radiation field for a half time step Δt/2 as described in the previous section to determine the energy and momentum source terms S_r (E) and S _r ( P ) (Equation (30)). Add these terms to the gas momentum and total energy.
• (d)  
  Apply the MHD solver for a full time step Δt using the partially updated variables at Δt/2.
• (e)  
  Repeat the implicit solver of the radiation field for a full time step Δt and add the energy, as well as momentum source terms, to the gas.
• (f)  
  Continue the above steps until the end of the simulation.

To test our treatment of the source term, we set up a uniform box in 2D covering the region [0, 1] × [0, 1]. The box is initialized with a uniform density ρ = 1 and gas temperature T = 1. The radiation field is initialized isotropically with the mean energy density E_r = 100. We fix the two dimensionless parameters ${\mathbb{P}}$ and ${\mathbb{C}}$ to be 1 and 100, respectively. We use 32 grid points for both directions. A constant absorption opacity ρ κ_a = 100 is adopted for the test. The typical thermalization timescale is then $1/\left({\mathbb{C}}\rho {\kappa }_{a}\right)={10}^{-4}$, while the time step we use is roughly 0.001. Therefore, in one time step, the system reaches the correct thermal equilibrium state and stays there, as shown in the left panel of Figure 1. The equilibrium state also agrees with the analytical solution very well as constrained by the total energy conservation. During the thermalization process, the algorithm can also keep the sign of E_r − a_r T⁴ so that there is no oscillatory behavior as a function of time, which can happen for high-order time integrators when the time step is much larger than the thermalization timescale. As another test, we change the initial condition to be T = 100, E_r = 1, and ρ κ_a = 1 so that the gas is much hotter than the equilibrium state in this case. We can still get the correct equilibrium state even though there is a huge difference between a_r T⁴ and E_r initially, as shown in the right panel of Figure 1.

Zoom In Zoom Out Reset image size

Our treatment of the source term is very similar to the scheme as described in Jiang et al. (2014), where the flux divergence term is calculated explicitly. However, there are nontrivial differences for the term n · ∇ I, even for the case with a uniform spatial distribution. When this term is calculated using the variables from the last time step, it is guaranteed to be zero to the level of the round-off error, as there is no spatial gradient at each step. This is not the case for the fully implicit scheme. During each cycle of the iteration when the gradient is calculated for each cell, variables in the neighboring cells are taken from the last iterative step while we use the values from this step for the current cell. Before the system reaches the equilibrium state, the flux divergence term is actually not zero. After the solution converges, we recover the same solution as given by the explicit scheme with the accuracy set by the tolerance level. The Jacobi-like iterative scheme we adopt also guarantees that all of the cells are treated in the same way during each iteration. Therefore, we do not generate any artificial spatial gradient during the iteration.

We have also tested the velocity-dependent terms by giving the gas a uniform horizontal velocity, v_x = 3. Specific intensities are initialized isotropically in the lab frame with the mean energy density E_r = 1. The initial gas temperature and density are T = 1 and ρ = 1, respectively. We choose $\rho {\kappa }_{a}=1,{\mathbb{P}}=1$ and ${\mathbb{C}}=100$ for this test and use 12 angles per cell. The system settles down to a steady state with gas velocity reduced to v_x = 2.956 as a small fraction of gas momentum is transferred to photons. The horizontal lab frame radiation flux in the steady-state F_r,x is very close to v_x (E_r + P_r,xx ), which is the value given by the first-order expansion in terms of v/c in radiation moment equations (Lowrie et al. 1999), since ${v}_{x}/{\mathbb{C}}$ is only 3% in this test. During this process, the change of kinetic energy is converted to radiation energy and lab frame E_r = 1.13 in the steady state. The specific intensities can also be compared with the analytical values directly, which are shown in Figure 2. In the steady state, the radiation field is isotropic in the comoving frame. Angular distribution of specific intensities in the lab frame is then determined by the Lorentz transformation (Equation (2)). The numerical solution agrees with the analytical values very well.

Zoom In Zoom Out Reset image size

Notice that the radiation field is not isotropic in the lab frame, and the Eddington tensor is P_r,xx /E_r = 0.417. This will be different from the Eddington tensor returned by the VET method for this case, since the short characteristic algorithm typically neglects the velocity-dependent terms, and it will return an Eddington tensor 1/3 I in the lab frame.

Our iterative scheme is typically easier to converge when the source term dominates over the transport term, because this means the diagonal coefficients in the resulting matrix of Equation (20) are larger than the off-diagonal one. The most difficult scenario is the transport of photons through a medium with zero opacity, where all of the source terms are zero and we only iterate over the transport term. To show how the algorithm works in this case, we inject two beams of photons from the bottom of a box covering the region [−0.5, 0.5] × [−2, 2]. At x = ±0.1, y = −2, the specific intensities with μ_x = ±0.577, μ_y = 0.577 are set to be 0.8, while the specific intensities along the other angles are zero. All of the intensities inside the box are also set to zero initially. The box is periodic along the x-direction and outflow at the top. The opacity coefficients of the gas are set to zero, and the gas does not evolve in this test. The spatial resolution is 64 × 256 for the horizontal and vertical directions, and we use four angles per cell. The time step is 4 × 10⁻³, which is also the light-crossing time across the vertical direction of the box. During the first step, it takes 100 iterations to reach the relative error of 10⁻³, which declines very slowly with more iterations. After one step, the two beams have crossed the whole simulation box. In about 10 time steps, the system reaches a steady state, and the spatial distribution of E_r is shown in the left panel of Figure 3. Compared with the same test done by Jiang et al. (2014), the solution is very similar to the case with first-order spatial reconstruction but more diffusive compared with the solution using second-order spatial reconstruction. This is expected, as our implicit algorithm only uses first-order spatial reconstruction. The overall cost of the implicit algorithm is also comparable to the cost if the transport equation is solved explicitly for this case, because the number of iterations we need for convergence is now comparable to the ratio between the speed of light and the sound speed. Another reason that this test requires a lot of iterations for convergence is that the initial condition is far away from the steady-state solution. The injected photons at the bottom boundary need to be communicated with the whole simulation box during the iterations. However, a multigrid technique can significantly speed up the convergence in this case because the information at the boundary can propagate to the domain much faster at lower resolution levels.

Zoom In Zoom Out Reset image size

Our algorithm is a finite volume scheme, which evolves the specific intensities using fluxes at cell faces as other cell-centered variables such as density in Athena++. Therefore, the existing mesh refinement module in the code (Stone et al. 2020) can be directly used for our RT algorithm. One example is shown in the right panel of Figure 3. The setup is the same as the left panel except that we refine the region [−0.5, 0.5] × [−1, 1] with two levels of refinements. Photons can cross the refinement boundaries smoothly without any numerical artifacts.

One widely used test to demonstrate the difference between solving the angular resolved RT and radiation moment equations based on simplified assumptions of the closure relation is to capture the shadow cast by an optically thick cloud with beams of photons (Hayes & Norman 2003; González et al. 2007). Diffusion-like approximation will always send photons to the region where the radiation energy density gradient exists and thus cannot keep the shadow. Other approximations that treat the radiation field as a fluid (such as M1; González et al. 2007; McKinney et al. 2014) have their own caveats, such as artificially merging photons coming from different directions. Our implicit scheme should get the correct shadows, as we still solve the specific intensities along different directions independently.

We used the same setup as in Jiang et al. (2012, 2014) for direct comparison. The simulation domain covers the Cartesian box [−0.5, 0.5] × [−0.3, 0.3] with resolution 512 × 256. We use 40 angles per cell for the 2D simulation. The gas density has the distribution $\rho (x,y)=1+9/\left[1+\exp \left[10\left({\left(x/0.1\right)}^{2}+{\left(y/0.06\right)}^{2}-1\right)\right]\right]$, while the temperature T is set to achieve a constant gas pressure P = 1 for the whole box. We use pure absorption opacity κ_a = ρ T^−3.5 so that the total optical depth across the box in the hot gas around the cloud is only 1, while the optical depth across the cloud is more than 10⁵. For the left boundary, specific intensities along the directions ±15° with respect to the horizontal axis have a fixed value of 103.13, but all others are set to zero. For the right boundary, outgoing specific intensities are just copied from the last activated zones to the ghost zones, while incoming specific intensities are zero. A periodic boundary condition is used for the top and bottom boundaries. We choose the dimensionless speed of light to be ${\mathbb{C}}=1.9\times {10}^{5}$, which does not affect the steady-state solution. All specific intensities are initialized to zero, and we keep the gas quantities fixed for the shadow test. It takes about 1000 iterations to reach a relative error of 10⁻⁶ with three time steps, and the radiation quantities are very close to the steady-state values across the whole simulation box. The spatial distributions of the mean radiation energy density E_r and the components of the Eddington tensor f_xx and f_yy are shown in Figure 4. Our implicit algorithm gets the umbra and penumbra correctly, and the results are very similar to the solutions returned by the explicit scheme as shown in Jiang et al. (2014), as well as the VET method as described by Jiang et al. (2012).

Zoom In Zoom Out Reset image size

One great advantage of our implicit algorithm is that it can also capture the time-dependent evolution of the radiation field efficiently, even though the speed of light is 1.9 × 10⁵ times the typical sound speed. To demonstrate that our algorithm can also capture the full dynamical evolution of the cloud, we use the same setup as in the shadow test but evolve the hydrodynamic quantities with the radiation field together. We choose the temperature scale such that ${\mathbb{P}}={10}^{-7}$ so that the radiation pressure of the injected photons from the left boundary is about 10⁻³ of the initial gas pressure. Notice that the static shadow test is independent of this ratio. The parameter is also chosen to match the same experiment done by Jiang et al. (2012). An outflow boundary condition is used for the hydrodynamic variables at the left and right boundaries, while a periodic boundary is used for the top and bottom. The typical sound-crossing timescale across the box is t₀ = 1, while the light-crossing time in the low-density gas is ≈ 5.2 × 10⁻⁶ t₀. The photon diffusion timescale across the cloud is comparable to t₀. Since the effective temperature of the incoming radiation field is larger than the temperature of the cold gas by a factor of 30, the gas will be heated up in the timescale of 2 × 10⁻⁵ t₀. The density distributions of the cloud at three different times t = 0.06t₀, 0.1t₀, and 0.16t₀ are shown in Figure 5. The cloud is initially heated up at the surface of the left side. The heated gas drives shocks into the cloud, and some gas is also evaporated simultaneously. The shock propagates inside the cloud and eventually gets destroyed. The overall evolution is very similar to the results shown by Jiang et al. (2012) and Proga et al. (2014), despite some small differences of the detailed structures. Compared with the VET scheme, there is a subtle difference in the equations that are solved. The VET scheme typically calculates the Eddington tensor at the beginning of each time step, which is then held fixed when the radiation moment equations and hydrodynamic equations are evolved for one time step. The algorithm developed here solves the specific intensities implicitly with the same time evolution as the hydrodynamic variables. If we integrate the specific intensities over the angular space, we effectively evolve the radiation moments using the Eddington tensor determined at the end of each time step, consistent with the implicit update of the radiation field. This will not make any difference for steady-state solutions or if the time step is sufficiently small. However, when the time step is much larger than the light-crossing time for each cell, a change of the Eddington tensor due to time evolution of the radiation field within one time step may not be negligible. Future studies are needed to quantify the difference.

Zoom In Zoom Out Reset image size

Photon diffusion in the optically thick regime, which is easy to solve with the classical diffusion approximation by design, turns out to require careful treatment when the full time-dependent transport equation is solved. This is because in the very optically thick regime, the numerical diffusion associated with the standard HLLE solver can be much larger than the physical diffusive flux (see the Appendix of Jiang et al. 2013). To demonstrate that our treatment of the transport term as described in Section 3.2.1 is sufficient to capture the photon diffusion correctly, we set up a radiation field with mean energy density ${E}_{r}(x)=\exp \left(-40{x}^{2}\right)$ in the box $x\in \left(-0.5,0.5\right)$. For ∣x∣ > 0.5, we set ${E}_{r}=\exp \left(-10\right)$. We use 256 grid points for the whole simulation box (−1, 1) and two angles per cell so that the Eddington tensor is always one-third. An outflow boundary condition is used for this test. The gas has a constant scattering opacity ρ κ_s = 4 × 10⁴, and we do not evolve the gas quantities for this test. We choose the dimensionless speed of light to be ${\mathbb{C}}=10$. The evolution of E_r can be described by the solution to the diffusion equation as (Sekora & Stone 2010; Jiang et al. 2014)

$$\begin{eqnarray}&&{E}_{r}(x,t)=\displaystyle \frac{1}{{\left(160{Dt}+1\right)}^{1/2}}\exp \left(\displaystyle \frac{-40{x}^{2}}{160{Dt}+1}\right),\end{eqnarray} \tag{ 33 }$$

where $D\equiv {\mathbb{C}}/(3\rho {\kappa }_{s})$ is the diffusion coefficient. The specific intensities are initialized isotropically with the mean energy density matching the initial profile of E_r . Profiles of E_r at three different times from the numerical solution are shown in Figure 6 and match the analytical solution very well. Notice that because we do not set the boundary condition to match the diffusion solution exactly, we can only compare the region roughly $\sqrt{{Dt}}$ away from the boundary. We have also tried cases with larger opacity, and our algorithm can still get the correct solution. However, we need to use a larger value of the parameter τ_c (Equation (17)) to get the same accuracy with a larger κ_s .

Zoom In Zoom Out Reset image size

Another important regime is dynamic diffusion, where the photons are advected with the flow while diffusion happens. We use a similar setup as for the static diffusion case but extend the simulation domain to (−10, 10) with 1280 grid points. A periodic boundary condition is used so that the photons are not advected out of the box. We give the flow a constant velocity, v = 1. The dimensionless speed of light is also increased to ${\mathbb{C}}=1000$ for this test so that the typical diffusion speed is not negligible compared with the advection speed. Profiles of the mean radiation energy densities at time t = 4, 8, and 16 are shown as the solid lines in Figure 7. The analytical solutions can also be estimated by replacing the position x in Equation (33) with x − vt to the first order of v/c, which is 10⁻³ in this test. The analytical solutions are shown as the black dashed lines in Figure 7, and they match the numerical solution very well.

Zoom In Zoom Out Reset image size

The test of a uniform temperature, scattering-dominated atmosphere done by Davis et al. (2012) is useful to show how our iterative scheme performs when there is a large dynamical range of radiation energy density with both optically thin and optically thick regions. In particular, we can check the performance of our iterative scheme for problems that can be solved efficiently with the classical short characteristic method. We use 1280 grid points covering a simulation domain x ∈ [−10, 10] in 1D. The gas density follows the exponential profile $\rho ={10}^{-3}\exp (10-x)$, and the gas temperature is fixed to be T_i = 1. The ratio between absorption and total opacity is , which is taken to be a constant in the whole simulation domain. Therefore, the scattering and absorption opacity are κ_a = and κ_s = 1 − . The gas quantities are fixed to be the initial conditions with zero velocity and do not evolve. Since our algorithm always solves the gas internal energy and radiation field together, we simply reset the gas quantities to be the initial profiles at the end of each step. We use one angle per octant for this test so that the Eddington tensor is always one-third. At the bottom boundary, we set the specific intensities to be the isotropic value ${a}_{r}{T}_{i}^{4}/(4\pi )$. At the top boundary, the incoming specific intensity is set to zero, while the outgoing specific intensity is copied from the last active zone to the ghost zone. Specific intensities inside the simulation domain are initialized to be ${a}_{r}{T}_{i}^{4}/(4\pi )$.

The steady-state profile of radiation energy density for this atmosphere setup is

$$\begin{eqnarray}&&{E}_{r}={a}_{r}{T}_{i}^{4}\left[1-\displaystyle \frac{\exp \left(-\sqrt{3\epsilon }\tau \right)}{1+\sqrt{\epsilon }}\right],\end{eqnarray} \tag{ 34 }$$

where $\tau \equiv {\int }_{x}^{10}\rho {dx}$ is the total optical depth integrated from x = 10. Upward and downward specific intensities are given by ${I}_{+}=(1/4\pi )\left[{E}_{r}+(1/3{\mu }_{x})\partial {E}_{r}/\partial \tau \right]$ and ${I}_{-}=(1/4\pi )\left[{E}_{r}-(1/3{\mu }_{x})\partial {E}_{r}/\partial \tau \right]$, respectively. The steady-state numerical solutions for varying from 0.1 to 10⁻⁸ are shown as the solid lines in Figure 8, which agree with the analytical solutions very well. Our algorithm requires more iterations to reach the same level of accuracy when decreases. The relative error is also larger in the optically thin surface, and the error increases with reduced for a fixed spatial resolution. All of these properties are very similar to the results given by the short characteristic method for this test problem (Davis et al. 2012). We have also tried the test in a 3D domain with the atmosphere profile aligned with the x-axis. We use a periodic boundary condition for the y- and z-directions in this case. The steady-state solutions are identical to the 1D case, and the rate of convergence is also the same. The overall cost is just proportional to the total number of grid points.

Zoom In Zoom Out Reset image size

To quantify the rate of convergence for our algorithm and compare with the short characteristic method, we pick the = 10⁻⁶ case and measure the relative error ΔI/I of the specific intensities at x = 10 after each iteration. Here we calculate ΔI as the difference between the solution at each iterative step and the analytical solution. We only focus on the location at x = 10 because the difference between the initial condition and the analytical solution is the largest there. It is also the most optically thin region and takes more iterations to converge. The criterion is different from what Davis et al. (2012) used because we do not iterate over the source term. But both criteria serve the same purpose. Notice that the steady-state solution does not depend on the speed of light. However, our algorithm always solves the time-dependent transport equation and relaxes to the steady-state solution with the time step Δt fixed by the gas sound speed v₀. We have tried three different values of c/v₀, and the changes of ΔI/I with the number of iterations are shown in Figure 9. When c/v₀ is small, although the implicit scheme is easier to converge per time step, it actually takes a larger total number of iterations to get the steady-state solution. When c/v₀ is larger than 10⁴ so that cΔt is larger than the length of the simulation box, the convergence rate is independent of c/v₀. Our implicit algorithm also takes more iterations to converge with this initial condition compared with the short characteristic method. This test demonstrates that the implicit scheme is not optimal just to find the steady-state solution, as we are solving the full time-dependent RT, although we can still get the correct answer.

Zoom In Zoom Out Reset image size

Solving the RT equations in 1D spherical polar coordinates can be useful for many systems that are close to spherically symmetric (such as stars), particularly with the future development of multiple frequency groups. Spherical polar angular systems as described in Section 3.2.4 need to be used in this case because of the assumed symmetry. One useful test for this case is the homogeneous sphere problem, which is widely used to test the numerical scheme for radiation and neutrino transport (Abdikamalov et al. 2012; Radice et al. 2013). For this test, we set up the gas with a constant density ρ₀ = 1 and temperature T₀ = 1 from radius R₀ = 0.05 to 1. Between R₀ = 1 and the outer boundary, which is chosen to be 7, the density is 10⁻⁷, and the temperature is 10⁻³. We assume an absorption opacity ρ₀ κ_a = 10 for R₀ smaller than 1, and the opacity is zero in other regions. The gas is held with zero velocity and does not evolve. The units can be chosen to match any system and do not affect the solution to the transport equation. We set the specific intensity to be isotropic, with the value ${a}_{r}{T}_{0}^{4}/(4\pi )$ at the inner boundary. For the outer boundary, we require the outgoing specific intensities to be continuous and the incoming specific intensities to be zero. We use 1000 uniformly distributed grid points for the spatial resolution and 40 angles per cell for the specific intensities. This setup covers both the optically thick and a sharp transition to the optically thin region, which can be challenging if the transport equation is not solved accurately. The specific intensities are initialized to be the isotropic value ${a}_{r}{T}_{0}^{4}/(4\pi )$ in the whole grid. We choose the speed of light to be 100 times the gas sound speed, and the code finds the steady-state solution in roughly 20 iterations. Notice that the steady-state solution does not depend on the speed of light.

This setup also has an analytical solution, which we can compare with. The specific intensities as a function of radius r and angles n are (Smit et al. 1997)

$$\begin{eqnarray}&&I(r,{\boldsymbol{n}})=\displaystyle \frac{{a}_{r}{T}_{0}^{4}}{4\pi }\left[1-\exp \left({\rho }_{0}{\kappa }_{a}s(r,{n}_{r})\right)\right],\end{eqnarray} \tag{ 35 }$$

where n_r is the radial component of the unit vector n and the function s(r, n_r ) is defined as

$$\begin{eqnarray}s(r,{n}_{r})=\left\{\begin{array}{ll}{{rn}}_{r}+{R}_{0}g(r,{n}_{r}) & \mathrm{if}\ r\lt {R}_{0}\\ 2{R}_{0}g(r,{n}_{r}) & \mathrm{if}\ r\geqslant \ {R}_{0}\ \&\\ & \sqrt{1-{\left({R}_{0}/r\right)}^{2}}\leqslant {n}_{r}\leqslant 1\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 36 }$$

with the function g(r, n_r ) to be

$$\begin{eqnarray}&&g(r,{n}_{r})=\sqrt{1-{\left(\displaystyle \frac{r}{{R}_{0}}\right)}^{2}\left(1-{n}_{r}^{2}\right)}.\end{eqnarray} \tag{ 37 }$$

Moments of the radiation field (E_r , F_r , P_r ) can be calculated directly via angular quadrature of I(r, n ).

Radial profiles of E_r , F_r , and P_r /E_r from the numerical solution are shown as the solid lines in Figure 10 and agree with the analytical solutions very well. In particular, our numerical algorithm can capture the sharp transition at r = 1 without showing any artificial numerical effect. The Eddington tensor P_r /E_r is one-third deep in the sphere. It is slightly smaller than one-third around r = 1. At a large distance, the sphere effectively becomes a point source, and only outgoing specific intensities pointing radially are nonzero. Therefore, the Eddington tensor approaches 1. The slight difference between the numerical and analytical solutions for P_r /E_r is due to the finite angular resolution to resolve the angular domain $\sqrt{1-{\left({R}_{0}/r\right)}^{2}}\leqslant {n}_{r}\leqslant 1$ at large r.

Zoom In Zoom Out Reset image size

Properties of convergence for the full radiation MHD algorithm can be checked with linear wave tests. The linearized equations we solve here are identical to what were evolved by Jiang et al. (2014). In order to compare directly, we adopt the same setup here. The background gas has a uniform density ρ = 1 and temperature T = 1 with adiabatic index γ = 5/3. We use one angle per octant with a mean radiation energy density E_r = 1 so that the Eddington tensor is always 1/3 I , as used in the analytical solutions. We fix the dimensionless speed of light to be ${\mathbb{C}}=10$ and only use a constant pure absorption opacity. We use a 2D simulation box with size L_x = 1, L_y = 1 and change the spatial resolution from 32 × 32 to 512 × 512. Periodic boundary conditions are used for both directions. We measure the L1 error between the numerical and analytical solutions over the whole simulation box for each resolution and check how the errors change with resolution. We have also tried the same test in 3D and get the same results, although these 3D tests are more expensive to run. We initialize the calculation with the eigenmodes of the radiation-modified sound waves as given by Appendix B of Jiang et al. (2012) and stop the calculation after the waves propagate one period along the x-direction.

The two key dimensionless parameters that control the properties of radiation-modified sound waves are optical depth per wavelength length τ_a and the ratio between radiation pressure and gas pressure ${\mathbb{P}}$ in the background state after ${\mathbb{C}}$ is fixed. Figure 11 shows the change of with resolution N for three different values of τ_a and two different values of ${\mathbb{P}}$. When the error is dominated by the hydrodynamic part, as in the optically thin regime, small radiation pressure, or low resolution, the error decreases with increasing resolution as N². When the error is dominated by the RT module, it only shows first-order convergence, as expected. Part of the truncation error in RT is due to different upwind directions for specific intensities propagating along opposite directions, as explained in Section 3.2.1. For the case with the highest resolution, N = 512, we find that increasing τ_c (Equation (17)) in the HLLE flux by a factor of 2 from our default value can help decrease the L1 error.

Zoom In Zoom Out Reset image size

We can also measure the numerically calculated propagation speed and damping rate of the linear waves with the resolution N = 256, which are shown in Figure 12 for τ_a varying from 10⁻² to 10² with three different values of ${\mathbb{P}}$ in each case. When the optical depth is small, our algorithm is able to capture both the adiabatic sound wave with ${\mathbb{P}}=0.01$ and the isothermal sound wave with ${\mathbb{P}}=10$. When τ_a is larger than 10 so that radiation and gas are tightly coupled, we are also able to capture the radiation-driven acoustic modes correctly. The results are the same as shown in Figure 10 of Jiang et al. (2014), as the same parameters are used. But the time step in this algorithm is a factor of 10 larger compared with the case when the time step is limited by the speed of light.

Zoom In Zoom Out Reset image size

Another challenging test for the full algorithm is the radiation-modified shocks, which check the numerical scheme in the nonlinear regime. In particular, we can see how well the source terms balance the transport terms for both the RT and hydrodynamic parts, which is crucial to get the structures of the shock correct, even though the time step is much larger than the light-crossing time per cell. Since our algorithm has the full angular distribution of the radiation field, we will compare with the semianalytical solutions provided by Ferguson et al. (2017), which gets the radiation shock solutions by solving the time-independent RT and hydrodynamic equations numerically without making any assumption on the closure relation. These are the same radiation shock solutions that are used to test the explicit RT algorithm in Jiang et al. (2014).

We choose a unit system such that the upstream adiabatic sound speed is 1 and the dimensionless speed of light ${\mathbb{C}}=1.73\times {10}^{3}$. The upstream gas has density and temperature ρ₀ = 1 and T₀ = 0.6 with adiabatic index γ = 5/3. We change the upstream flow velocity for cases with different Mach numbers. The dimensionless parameter ${\mathbb{P}}$ is 7.72 × 10⁻⁴, which is roughly the ratio between the radiation pressure and gas pressure in the upstream. We choose a constant absorption opacity ρ κ_a = 577.4 for the whole simulation domain. We initialize the 1D simulation domain with the semianalytical solutions given by Ferguson et al. (2017) and see how well the code can keep the shock structures. We use 4096 grid points for all of the tests. For the left boundary, all of the gas quantities are fixed to be the upstream values, while the specific intensities are isotropic in the comoving frame with a mean energy density fixed to be ${a}_{r}{T}_{0}^{4}$. For the right boundary, we simply copy the gas quantities from the last active zone to the ghost zones, but we find that it is necessary to keep the gradients of specific intensities continuous across the boundary.

Profiles of various quantities after the numerical solution reaches a steady state are shown in Figure 13 for the case with upstream flow velocity v = 2. This is the case with a subcritical shock (Mihalas & Mihalas 1984; Sincell et al. 1999; Lowrie & Edwards 2008), where the gas temperature has a Zel'dovich spike (Zel'Dovich & Raizer 1967) while the radiation temperature is continuous. For the subcritical case, the gas temperature before the spike is smaller than the downstream temperature. Particularly, the Eddington tensor is smaller than one-third on the right-hand side of the spike and larger than one-third on the left-hand side of the spike. Our numerical solution in the steady state is identical to the semianalytical solution, shown as the dashed red lines. Notice that our time step is set by the flow speed, and we only need roughly 20 iterations per time step. The numerical cost is significantly smaller than the explicit scheme in Jiang et al. (2014) for the same test, and it is also more efficient than the VET scheme for a similar test done by Jiang et al. (2012). The steady-state numerical solution with upstream flow velocity increased to v = 3 is shown in Figure 14. This is the case with a supercritical shock, which means the gas temperature at the upstream side of the spike is the same as the downstream temperature. The downstream gas is hotter, with an increased ratio between radiation pressure and gas pressure. Our algorithm is still able to keep the shock structures very well, as demonstrated in Figure 14, by comparing our numerical solution with the semianalytical solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulation of radiation-driven dusty wind mentioned in the Introduction is a useful test to demonstrate the difference between different approaches to solve the RT equation. This is a test that requires solving the full radiation hydrodynamic equations, and the simulation results can be significantly different if the RT equation is not solved accurately. In order to compare with the results based on the VET method directly, we use the same setup as in the run T3_F0.5 done by Davis et al. (2014). The simulation box is a 2D Cartesian domain covering the region [−256, 256]L₀ × [0, 1024]L₀, where the length scale L₀ = 2.01 × 10¹⁵ cm is the fiducial scale height. The spatial resolution is 512 × 1024, and we use 40 angles per cell. The fiducial temperature is chosen to be T₀ = 82 K, with a corresponding fiducial velocity v₀ = 5.4 × 10⁴ cm s⁻¹ and fiducial timescale t₀ = 3.72 × 10¹⁰ s. We adopt the same Rosseland and Planck mean dust opacities as in Davis et al. (2014):

$$\begin{eqnarray}\begin{array}{rcl}{\kappa }_{a} & = & 0.0316{\left(\displaystyle \frac{T}{10{\rm{K}}}\right)}^{2}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},\\ {\kappa }_{{aP}} & = & 0.1{\left(\displaystyle \frac{T}{10{\rm{K}}}\right)}^{2}\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}.\end{array}\end{eqnarray} \tag{ 38 }$$

Notice that κ_{δ P} = κ_aP − κ_a is what we need in Equation (6). We initialize the gas with a constant density ρ₀ = 7.02 × 10⁻¹⁶ g cm⁻³ for 9 < y/L₀ < 10 and density floor 10⁻⁸ ρ₀ for all other regions. The gas has a uniform temperature T = T₀ and a constant gravitational acceleration g = 1.45 × 10⁻⁶ cm s⁻² along the −y direction. The radiation field is initialized to give ${E}_{r}={a}_{r}{T}_{0}^{4}$ and F_r = 0.5c² g/κ_a (T₀) through the whole simulation box. We assume that upward specific intensities have the same value at each cell, and the same assumption is made for downward specific intensities. The same approach is used to set the bottom boundary condition with the same F_r but ${E}_{r}=9.26{a}_{r}{T}_{0}^{4}$. This is designed to maintain a constant vertical radiation flux F_r with an initial total optical depth τ₀ = 3. At the top boundary, we set incoming specific intensities to zero and just copy the outgoing specific intensities from the last active zones to the ghost zones. For gas quantities, we use a reflecting boundary condition at the bottom and an outflow boundary condition at the top. We use a random perturbation on density with amplitude 12.5% to seed the instability. The ratio between the speed of light and the fiducial sound speed is ${\mathbb{C}}=5.54\times {10}^{5}$ for this test.

Density distribution at two snapshots t = 35t₀ and 80t₀ are shown in Figure 15. The gas shell is initially pushed away by the radiation force, and Rayleigh–Taylor instability has developed at t = 35t₀. The high-density filaments fall back temporarily, but they get accelerated again later on. At time t = 80t₀, the majority of the gas is still flowing outward. At each snapshot, we calculate the box-averaged Eddington ratio f_E , volume (τ_V ), and flux-weighted (τ_F ) optical depths as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{E} & = & \displaystyle \frac{\langle \rho {\kappa }_{a}{F}_{r}\rangle }{{cg}\langle \rho \rangle },\\ {\tau }_{V} & = & 1000{L}_{0}\langle {\kappa }_{a}\rho \rangle ,\\ {\tau }_{F} & = & 1000{L}_{0}\displaystyle \frac{\langle {\kappa }_{a}\rho {F}_{r,y}\rangle }{\langle {F}_{r,y}\rangle },\end{array}\end{eqnarray} \tag{ 39 }$$

where 〈 · 〉 represents the volume average over the whole simulation box. Histories of f_E , τ_V , and τ_F /τ_V are shown in Figure 16. The Eddington ratio dips below 1 around 30t₀ due to the Rayleigh–Taylor instability but returns to 1 quickly. The timescale when the dip shows up and the evolution history are basically the same as in Figure 6 of Davis et al. (2014), despite some small difference at the early phase, which is likely due to slightly different initial conditions. This demonstrates that the two completely different algorithms capture the same growth rate and nonlinear structure of the radiation-driven Rayleigh–Taylor instability. The Eddington ratio in our simulation eventually becomes larger than 1, and the gas is blown away from the top boundary. The 2D density structures shown in Figure 15 can also be directly compared with Figure 4 in Davis et al. (2014), and they share very similar structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is also interesting to compare the performance of our algorithm for this test with the performance of the VET scheme. The time step for our algorithm is 6 × 10⁻³ t₀, and we use 50 iterations per time step to achieve a relative accuracy of 10⁻⁶, except for the first few time steps, where 200 iterations are used per time step to pass the initial transient. With 256 MPI ranks using Skylake nodes, we are able to run roughly 2 × 10⁴ steps to reach 120t₀ within about 4 hr of wall-clock time. The cost per time step is comparable to the cost of the VET method used by Davis et al. (2014). However, Davis et al. (2014) used a time step that was a factor of 10 smaller to enable convergence, which made the overall cost significantly larger than the cost of our new scheme developed here. Our implicit scheme also reduces the cost by a factor of ≈10⁴ if we solve the problem using an explicit scheme by limiting the time step with the speed of light.