AN ALGORITHM FOR RADIATION MAGNETOHYDRODYNAMICS BASED ON SOLVING THE TIME-DEPENDENT TRANSFER EQUATION

•   
  Step 1. Calculate time step Δ t based on cell size and speed of light with Courant–Friedrichs–Lewy number 0.4.
•   
  Step 2. Calculate the change of specific intensity along each direction ${{\boldsymbol n}}_l$ due to divergence of transport flux along three directions ΔF_l, according to Section 4.2. Depending on whether the scattering opacity is larger or smaller than the absorption opacity in each cell, ΔF_l is added to the right-hand side of the matrix either in Step 5 or Step 4. This is necessary to keep the balance between the transport term and source terms.
•   
  Step 3. Estimate flow velocity for half time step according to Section 4.5.
•   
  Step 4. Update all the specific intensities with the absorption opacity related terms according to Section 4.3. Calculate the corresponding radiation energy and momentum source terms.
•   
  Step 5. Repeat Step 4 but for the scattering opacity related terms.
•   
  Step 6. Update gas quantities using the usual Athena algorithm described in Stone et al. (2008).
•   
  Step 7. Add the radiation energy and momentum source terms calculated in Step 4 and Step 5 to the gas quantities.

The left-hand side of Equation (6) represents the advection of specific intensity at the speed of light along the direction specified by ${{\boldsymbol n}}$ when the source term vanishes. Based on this idea, Stone & Mihalas (1992) applied upwind monotonic interpolation methods to solve the transport step in one dimension. In this way, each I is advected across the grid using fluxes calculated at the cell interface, which guarantees conservation of radiation energy and momentum to roundoff error. A similar technique can be used in multi-dimensions with appropriate extensions.

In the non-relativistic case, the direction of the specific intensity can be assumed to be constant during the transport step. Therefore, the transport part of Equation (6) can be rewritten in conservative form

$$\begin{eqnarray} \frac{\partial I}{\partial t}+{\mathbb {C}}{{\boldsymbol \nabla}}\cdot ({{\boldsymbol n}}I) &=& \frac{\partial I}{\partial t}+{\mathbb {C}}\frac{\partial }{\partial x}(\mu _x I)\nonumber \\ &&+{\mathbb {C}}\frac{\partial }{\partial y}(\mu _y I)+{\mathbb {C}}\frac{\partial }{\partial z}(\mu _z I)=0.\qquad \end{eqnarray} \tag{ 11 }$$

This equation suggests that we can decompose the transport of I along direction ${{\boldsymbol n}}$ into three transport steps along each axis and each of them can be calculated with similar method as proposed by Stone & Mihalas (1992).

The original algorithm developed by Stone & Mihalas (1992) needs to be modified for the velocity dependent source terms in optically thick regime to reduce numerical diffusion, especially for the case with pure scattering opacity. In the Eulerian frame, the transport term ${\mathbb {C}}{{\boldsymbol \nabla}}\cdot ({{\boldsymbol n}}I)$ includes both the diffusive part and advective part caused by the flow velocity. The diffusion and advection speed can be much smaller than the speed of light. Therefore, if we always use the speed of light as the transport speed as originally adopted by Stone & Mihalas (1992), the amount of numerical diffusion can significantly overwhelm the physical flux. A similar issue also exists when the two radiation moment equations are solved in the mixed frame as discussed in the Appendix of Jiang et al. (2013b). Therefore, it is necessary to separate the diffusion and advection parts so that we can treat them accurately.

The velocity dependent source term that needs special treatment is

$$\begin{eqnarray} &&IV\equiv 3{{\boldsymbol n}}\cdot {{\boldsymbol v}}J. \end{eqnarray} \tag{ 12 }$$

We split the transport term as

$$\begin{eqnarray} &&{\mathbb {C}}{{\boldsymbol n}}\cdot {{\boldsymbol \nabla}}I={\mathbb {C}}{{\boldsymbol n}}\cdot {{\boldsymbol \nabla}}\tilde{I}+{{\boldsymbol n}}\cdot {{\boldsymbol \nabla}}(IV), \end{eqnarray} \tag{ 13 }$$

where $\tilde{I}\equiv I-IV/{\mathbb {C}}$. Mathematically, the equation is unchanged. However, the numerical treatments of the two terms are different to significantly reduce numerical diffusion. The transport velocity for the first term is chosen to be $\alpha {\mathbb {C}}$ (Jiang et al. 2013b) instead of ${\mathbb {C}}$, where

$$\begin{eqnarray} \alpha &=&\sqrt{\frac{1-\exp {(-\tau)}}{\tau }}, \nonumber \\ \tau &\equiv &[10\Delta l \times (\sigma _a+\sigma _s)]^2. \end{eqnarray} \tag{ 14 }$$

Here Δl is the cell size. When optical depth per cell is larger than 1, α ∼ 1/[10Δl(σ_a + σ_s)] and the transport speed is reduced by a factor of α. When the cell is optically thin, α ∼ 1 and the transport speed is unchanged. As discussed in Jiang et al. (2013b), the numerical factor 10 in the above equation is chosen so that the numerical diffusion does not affect the physical solution over a wide range of optical depth per cell in our tests, while still keeping the algorithm stable.

The speed of light is no longer the characteristic speed for the second term ${{\boldsymbol n}}\cdot {{\boldsymbol \nabla}}(IV)$ in Equation (13). This term represents the advection of specific intensity along direction ${{\boldsymbol n}}$ due to motion of the flow. Therefore, the transport speed for this term is flow velocity ${{\boldsymbol v}}$.

The two transport terms can both be decomposed along the three axes independently as

$$\begin{eqnarray} &&\frac{\partial I}{\partial t}+{\mathbb {C}}\frac{\partial }{\partial x}(\mu _x\tilde{I})+{\mathbb {C}}\frac{\partial }{\partial y}(\mu _y\tilde{I})+{\mathbb {C}}\frac{\partial }{\partial z}(\mu _z\tilde{I})=0, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \frac{\partial I}{\partial t} +\frac{\partial }{\partial x}(\mu _x IV)+\frac{\partial }{\partial y}(\mu _y IV)+\frac{\partial }{\partial z}(\mu _z IV)&=& 0. \end{eqnarray} \tag{ 16 }$$

Although the above two equations are no longer the simple advection equation solved by Stone & Mihalas (1992), similar techniques can still be used to calculate the transport flux, as they represent the same advection process with advection velocity $\alpha {\mathbb {C}}$ and $|{{\boldsymbol v}}|$ respectively, which are assumed to be constants during this step. The second-order van Leer scheme as Equation (4) of Stone & Mihalas (1992) is adopted to calculate the flux.

The absorption opacity related terms change the specific intensity according to the following equation

$$\begin{eqnarray} \frac{\partial I}{\partial t}&=&{\mathbb {C}}\sigma _a\left(\frac{T^4}{4\pi }-I\right) +3{{\boldsymbol n}}\cdot {{\boldsymbol v}}\sigma _a\frac{T^4}{4\pi } +{{\boldsymbol n}}\cdot {{\boldsymbol v}}\sigma _a I \nonumber \\ &&-\sigma _a\frac{{{\boldsymbol v}}\cdot {{\boldsymbol v}}}{{\mathbb {C}}}J- \sigma _a\frac{{{\boldsymbol v}}\cdot ({{\boldsymbol v}}\cdot {\sf K})}{{\mathbb {C}}}. \end{eqnarray} \tag{ 17 }$$

The zeroth and first angular moments of the above equation multiplied by 4π give the following equations for the change of radiation energy density and flux as

$$\begin{eqnarray} \frac{\partial E_r}{\partial t}&=&{\mathbb {C}}\sigma _a(T^4-E_r)+\sigma _a{{\boldsymbol v}}\cdot \left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right),\nonumber \\ \frac{\partial {{\boldsymbol F}}_r}{\partial t}&=&{-}{\mathbb {C}}\sigma _a\left({{\boldsymbol F}}_r-\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right)+{{\boldsymbol v}}\sigma _a(T^4-E_r).\qquad \end{eqnarray} \tag{ 18 }$$

The corresponding changes of gas temperature and momentum are

$$\begin{eqnarray} \frac{\rho R_{{\rm ideal}}}{\gamma -1}\frac{\partial T}{\partial t}&=&{-}{\mathbb {P}}{\mathbb {C}}\left(1-\frac{v^2}{{\mathbb {C}}^2}\right)\sigma _a(T^4-E_r) \nonumber \\ &&-2{\mathbb {P}}\sigma _a{{\boldsymbol v}}\cdot \left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right), \nonumber \\ \frac{\partial \rho {{\boldsymbol v}}}{\partial t}&=&{\mathbb {P}}\sigma _a\left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right)\nonumber \\ &&-{\mathbb {P}}\sigma _a(T^4-E_r)\frac{{{\boldsymbol v}}}{{\mathbb {C}}}. \end{eqnarray} \tag{ 19 }$$

Notice that the change of kinetic energy density due to the work done by radiation force is

$$\begin{eqnarray} \frac{\partial (\rho v^2/2)}{\partial t}&=&{\mathbb {P}}\sigma _a{{\boldsymbol v}}\cdot \left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right)\nonumber \\ &&-{\mathbb {P}}\sigma _a(T^4-E_r)\frac{v^2}{{\mathbb {C}}}. \end{eqnarray} \tag{ 20 }$$

Therefore, the change of total energy $\partial [{\mathbb {P}}E_r+\rho R_{{\rm ideal}}T/(\gamma -1)+\rho v^2/2]/\partial t=0$ as the source terms we add to the gas and radiation have exactly the same value but opposite signs. As the thermalization time scale ${\sim }1/({\mathbb {P}}{\mathbb {C}}\sigma _a)$ can be very short even compared with light crossing time per cell when the optical depth per cell is much larger than 1, Equation (17) must be solved implicitly to make the algorithm stable. When the flow velocity is not too small compared with speed of light ⁶ and radiation energy density is much larger than the gas pressure, the radiation work term is not negligible and should be solved simultaneously with other terms to get the correct equilibrium state. To avoid extra non-linear terms caused by the velocity-dependent terms, flow velocity is fixed during this step (see Section 4.5). The absorption opacity is also fixed during this step. To get the correct equilibrium state, change of the gas temperature also needs to be calculated simultaneously with Equation (17). Therefore, the specific intensity along each direction ${{\boldsymbol n}}_m$ ($I^{n+1}_m$) and gas temperature (T^{n + 1}) at time step n + 1 can be calculated based on their values at current time step n and time step Δt as

$$\begin{eqnarray} \frac{I_m^{n+1}-I_m^n}{\Delta t}&=&{\mathbb {C}}\sigma _a\left(\frac{(T^{n+1})^4}{4\pi }-I_m^{n+1}\right)\nonumber \\ &&+3{{\boldsymbol n}}_m\cdot {{\boldsymbol v}}\sigma _a\frac{(T^{n+1})^4}{4\pi } +{{\boldsymbol n}}_m\cdot {{\boldsymbol v}}\sigma _a I_m^{n+1} \nonumber \\ &&-\sigma _a\frac{v^2}{{\mathbb {C}}}\sum _l\left(I^{n+1}_lW_l\right) \nonumber \\ &&-\frac{\sigma _a}{{\mathbb {C}}}\sum _l\left\lbrace {{\boldsymbol v}}\cdot ({{\boldsymbol v}}\cdot {{\boldsymbol n}}_l{{\boldsymbol n}}_l)I^{n+1}_lW_l\right\rbrace,\nonumber \\ \frac{\rho R_{{\rm ideal}}}{\gamma -1}\frac{T^{n+1}-T^n}{\Delta t}&=&{-}8\pi {\mathbb {P}}\sigma _a{{\boldsymbol v}}\cdot \left\lbrace \sum _l \left({{\boldsymbol n}}_l I_l^{n+1}W_l\right)\right. \nonumber \\ &&\left. - \frac{1}{{\mathbb {C}}}\sum _l \left({{\boldsymbol v}}I^{n+1}_l W_l +{{\boldsymbol v}}\cdot {{\boldsymbol n}}_l{{\boldsymbol n}}_l I^{n+1}_lW_l\right)\right\rbrace \nonumber \\ &&-{\mathbb {P}}{\mathbb {C}}\left(1-\frac{v^2}{{\mathbb {C}}^2}\right)\sigma _a \left\lbrace (T^{n+1})^4 \right.\nonumber \\ &&\left. - 4\pi \sum _l\left(I_l^{n+1}W_l\right)\right\rbrace. \end{eqnarray} \tag{ 21 }$$

Here ranges for the lower scripts m and l are from 1 to the total number of angles in each cell. For total N angles for each cell, Equation (21) represents N + 1 equations, which we have to solve simultaneously. This set of equations are non-linear because of the (T^{n + 1})⁴ term. They are linear with respect to $I^{n+1}_m$. We use Newton–Raphson iteration to solve this set of equations with the initial guess chosen to be $I^n_m$ and Tⁿ. Because the non-linear terms are only fourth order polynomial, we usually find the iteration converges very quickly.

Newton–Raphson iteration requires inversion of the Jacobi matrix during each step, which can be simplified analytically to significantly reduce the computational cost. We first subtract the line of $I^{n+1}_N$ from all the lines for $I^{n+1}_m (m\ne N)$. Each $I^{n+1}_m (m\ne N)$ can now be expressed as a function of $I^{n+1}_N$ and T^{n + 1} easily. Therefore we only need to invert a 2 × 2 matrix for $I^{n+1}_N$ and T^{n + 1} numerically, after which all the $I^{n+1}_m (m\ne N)$ can be calculated directly. The total cost to invert the Jacobi matrix is only $\mathcal {O}(N)$.

With the updated solution $I^{n+1}_m$ and T^{n + 1}, the changes of the gas internal energy density ΔE_T and gas momentum $\Delta (\rho {{\boldsymbol v}})$ due to the absorption opacity terms are

$$\begin{eqnarray} &&\Delta E_T=\frac{\rho R_{{\rm ideal}}}{\gamma }(T^{n+1}-T^n), \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&\Delta (\rho {{\boldsymbol v}})=-\frac{{\mathbb {P}}}{{\mathbb {C}}}\left[\sum _l\left({{\boldsymbol n}}_l I_l^{n+1}W_l\right)- \sum _l\left({{\boldsymbol n}}_l I_l^{n}W_l \right) \right]. \end{eqnarray} \tag{ 23 }$$

The related change of kinetic energy density is

$$\begin{eqnarray} &&\Delta E_{K}=((\rho {{\boldsymbol v}}^n+\Delta (\rho {{\boldsymbol v}}))^2-(\rho {{\boldsymbol v}}^n)^2)/(2\rho). \end{eqnarray} \tag{ 24 }$$

We calculate the gas momentum source term $\Delta (\rho {{\boldsymbol v}})$ according to the conservation of gas and radiation momentum at time step n and n + 1, instead of calculating the momentum source term directly according to Equation (18). This not only guarantees momentum conservation to roundoff error but also avoids difficulties when absorption opacity σ_a is very large while ${{\boldsymbol F}}_r$ is very small. However, the sum of ΔE_T and ΔE_K differs from the change of radiation energy density on the order of $\mathcal {O}(\Delta t^2)$. In principle, this energy error can be reduced when we go to second-order accuracy in time, which requires solving Equation (21) twice. Our approach to calculate ΔE_T and ΔE_K assures that gas temperature is evolved exactly as Equation (21) describes and it will not be affected by the truncation error of the kinetic energy density. As we will show in the tests, the energy error is usually smaller than 0.01% and it will not accumulate.

The scattering opacity does not change the gas temperature. It only changes the flow velocity and thus the kinetic energy density. The equations we need to solve for the scattering opacity related terms are

$$\begin{eqnarray} \frac{\partial I}{\partial t}&=&{\mathbb {C}}\sigma _s(J-I) +{{\boldsymbol n}}\cdot {{\boldsymbol v}}\sigma _s(I+3J) -2\sigma _s{{\boldsymbol v}}\cdot {{\boldsymbol H}}\nonumber \\ &&+\sigma _s\frac{{{\boldsymbol v}}\cdot {{\boldsymbol v}}}{{\mathbb {C}}}J +\sigma _s\frac{{{\boldsymbol v}}\cdot ({{\boldsymbol v}}\cdot {\sf K})}{{\mathbb {C}}}. \end{eqnarray} \tag{ 25 }$$

Taking the zeroth and first moments of above equation and multiplying by 4π, we get the equations describing the evolution of radiation energy density and momentum as

$$\begin{eqnarray} \frac{\partial E_r}{\partial t}&=&{-}\sigma _s{{\boldsymbol v}}\cdot \left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right), \nonumber \\ \frac{\partial {{\boldsymbol F}}_r}{\partial t}&=&{-}{\mathbb {C}}\sigma _s\left({{\boldsymbol F}}_r-\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right). \end{eqnarray} \tag{ 26 }$$

The change of gas momentum due to scattering opacity is

$$\begin{eqnarray} &&\frac{\partial \rho {{\boldsymbol v}}}{\partial t}={\mathbb {P}}\sigma _s\left({{\boldsymbol F}}_r -\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}\cdot {\sf P_r}}{{\mathbb {C}}}\right). \end{eqnarray} \tag{ 27 }$$

The velocity-dependent terms in Equation (25) also needs to be calculated simultaneously with the other term in Equation (25) in order to get the correct equilibrium state, especially when the velocity dependent terms become significant. In order to make the code stable for the regime when the scattering optical depth per cell is larger than 1, the scattering opacity source terms also need to be added implicitly.

Following the approach used in Section 4.3 for the absorption opacity case, we fix flow velocity during this step to avoid the nonlinear terms. Then each specific intensity I_m is updated using backward Euler implicitly as

$$\begin{eqnarray} \frac{I_m^{n+1}-I_m^n}{\Delta t}&=&{\mathbb {C}}\sigma _s\left(\sum _l \left(I_l^{n+1}W_l\right)-I_m^{n+1} \right)\nonumber \\ &&+{{\boldsymbol n}}_m\cdot {{\boldsymbol v}}\sigma _s\left(I_m^{n+1}+3\sum _l\left(I_l^{n+1}W_l\right)\right)\nonumber \\ &&-2\sigma _s\sum _l\left({{\boldsymbol v}}\cdot {{\boldsymbol n}}_l I_l^{n+1}W_l\right) +\sigma _s\frac{v^2}{{\mathbb {C}}}\sum _l\left(I^{n+1}_lW_l\right)\nonumber \\ &&+\frac{\sigma _s}{{\mathbb {C}}}\sum _l\left({{\boldsymbol v}}\cdot \left({{\boldsymbol v}}\cdot {{\boldsymbol n}}_l{{\boldsymbol n}}_l I^{n+1}_lW_l\right)\right). \end{eqnarray} \tag{ 28 }$$

These N equations are linear with respect to the unknowns $I^{n+1}_m$ for given $I^{n}_m$, velocity and opacity. We find it is actually easier to use $I^{n+1}_mW_m$ as unknowns and the N × N matrix we need to invert has the following format

$$\begin{eqnarray} {\sf M_s} = \left[{\begin{array}{@{}ccccc@{}} a_1+b_1+c_1&b_2+c_1&\cdots & \cdots & b_N+c_1 \\ b_1+c_2&a_2+b_2+c_2&\cdots & \cdots & b_N+c_2 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ b_1+c_N & b_2+c_N & \cdots & \cdots & a_N+b_N+c_N \end{array}} \right]. \nonumber\\ \end{eqnarray} \tag{ 29 }$$

The matrix elements a_m, b_m and c_m are

$$\begin{eqnarray} a_m&=&\frac{1}{W_m}(1+\Delta t{\mathbb {C}}\sigma _s-\Delta t\sigma _s{{\boldsymbol n}}_m\cdot {{\boldsymbol v}}),\nonumber \\ b_m&=&{-}\Delta t{\mathbb {C}}\sigma _s+2\Delta t\sigma _s{{\boldsymbol v}}\cdot {{\boldsymbol n}}_m \nonumber \\ &&-\Delta t \sigma _s\frac{v^2}{{\mathbb {C}}}-\Delta t\frac{\sigma _s}{{\mathbb {C}}}({{\boldsymbol v}}\cdot ({{\boldsymbol v}}\cdot {{\boldsymbol n}}_l{{\boldsymbol n}}_l)),\nonumber \\ c_m&=&{-}3\Delta t{{\boldsymbol n}}_m\cdot {{\boldsymbol v}}\sigma _s. \end{eqnarray} \tag{ 30 }$$

LU decomposition is used to invert the matrix. Because of the special pattern of the matrix, the LU decomposition can be worked out analytically. The L and U matrix after the decomposition only have two lines that are non-zero. Therefore, the final cost to invert this matrix is still only $\mathcal {O}(N)$. After we get the updated solution $I^{n+1}_m$, we calculate the change of gas momentum and kinetic energy density according to Equations (23) and (24).

When source terms are added as described in Sections 4.3 and 4.4, the flow velocity is fixed to avoid the highly non-linear terms during the matrix inversions. When photon momentum is significant compared with gas momentum, this is no longer a good approximation. To improve the accuracy, instead of using the flow velocity at time step n, we first estimate the flow velocity at time step n + 1/2 based on the following equations:

$$\begin{eqnarray} \frac{\partial \rho {{\boldsymbol v}}}{\partial t}&=&{\mathbb {P}}(\sigma _s+\sigma _a)\left({{\boldsymbol F}}_r-\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}P_r}{{\mathbb {C}}}\right),\nonumber \\ \frac{\partial {{\boldsymbol F}}_r}{\partial t}&=&{-}{\mathbb {C}}(\sigma _s+\sigma _a)\left({{\boldsymbol F}}_r-\frac{{{\boldsymbol v}}E_r+{{\boldsymbol v}}P_r}{{\mathbb {C}}}\right), \end{eqnarray} \tag{ 31 }$$

where P_r is the diagonal components of the radiation pressure tensor ${\sf P_r}$. We keep E_r and P_r fixed in this step and solve ${{\boldsymbol v}}$ and ${{\boldsymbol F}}_r$ from the above two equations. For given values of ${{\boldsymbol v}}^{n}$ and ${{\boldsymbol F}}_r^{n}$, we estimate the velocity $\tilde{{{\boldsymbol v}}}$ according to

$$\begin{eqnarray} \rho \tilde{{{\boldsymbol v}}}-\rho {{\boldsymbol v}}^{n}&=&0.5\Delta t{\mathbb {P}}(\sigma _s+\sigma _a)\nonumber \\ &&\times\left(\frac{{\mathbb {C}}}{{\mathbb {P}}}\rho {{\boldsymbol v}}^n +{{\boldsymbol F}}_r^{n} -\frac{{\mathbb {C}}}{{\mathbb {P}}}\rho \tilde{{{\boldsymbol v}}}-\frac{\tilde{{{\boldsymbol v}}}}{{\mathbb {C}}}(E_r+P_r) \right).\quad \end{eqnarray} \tag{ 32 }$$

The estimated flow velocity $\tilde{{{\boldsymbol v}}}$ is used when the absorption and scattering source terms are added.

To demonstrate that our treatment of absorption opacity as described in Section 4.3 can give the correct solution for radiation energy density E_r and gas temperature T, we first study the evolution of E_r and T toward an equilibrium state in a static medium. It is important to show that even when the time step is much larger than the thermalization time, the method is still stable and accurate.

We setup a uniform 2D domain, periodic in x and y, with L_x = L_y = 1 and N_x = N_y = 32. Density is chosen to be 1 and flow velocity is zero. The dimensionless parameters are ${\mathbb {P}}=1, {\mathbb {C}}=10, R_{{\rm ideal}}=1$, and γ = 5/3. The typical sound speed is chosen to be only 10% of speed of light in order to speed up the calculation. We choose initial conditions with E_r ≠ T⁴ and follow it as it relaxes to a state with E_r = T⁴, subject to energy conservation. Figure 1 shows initial conditions E_r = 100, T⁴ = 1 for the left panel and T⁴ = 10⁸, E_r = 1 for the right panel. Evolution of E_r and T for the two different initial conditions with two different opacities σ_a are shown in Figure 1. The equilibrium solutions from the code agree with the analytical solutions (the blue lines) very well. The total energy error in this test is smaller than 10⁻¹⁰, which is set by the accuracy of the Newton–Raphson iterations as discussed in Section 4.3.

Zoom In Zoom Out Reset image size

When σ_a = 100, the time step, which is set by the speed of light, is much larger than the thermalization time. In this case, the correct equilibrium state is achieved within about one time step. When the time step is comparable to the thermalization time as in the case σ_a = 1, our algorithm automatically resolves the thermalization process. Unlike the predict and correct scheme shown in Figures 1 and 2 of JSD12, this algorithm can always keep the relative values of E_r and T⁴ even for extreme parameters during the thermalization process, as our algorithm solves the gas temperature and radiation energy density simultaneously.

To test our treatment of the velocity-dependent terms, we give the fluid a uniform initial velocity v_x along the x direction based on the setup used in the last section and let the system evolve. As the RT equations are not Galilean invariant (Lowrie et al. 1999, JSD12), specific intensity will be beamed along the direction of flow velocity and a net radiation flux is produced along the opposite direction. As the co-moving flux is non-zero initially, the flow experiences a drag by the photons and is decelerated until the co-moving flux becomes zero in the equilibrium state. We use $v_x=0.3{\mathbb {C}}$ in this test so that the velocity-dependent terms are significant. For this test, we set T = 1 and E_r = 1 initially. Three angles for each octant are used.

Histories of the flux, momentum, and energy for the case with absorption opacity σ_a = 100 are shown in Figure 2. The advective flux is defined to be $F_v\equiv v_x(E_r+ P_{r,xx})/{\mathbb {C}}$. When the system reaches the equilibrium state, F_v = F_r and the co-moving flux becomes zero, as it should be. The total momentum of the system is the sum of the gas momentum ρv_x and photon momentum ${\mathbb {P}}F_r/{\mathbb {C}}$, which is conserved as shown in the middle panel of Figure 2. The gas momentum is decreased due to the radiation drag and converted to photon momentum. During this process, the kinetic energy of the gas is transformed to the radiation energy density and gas internal energy through the thermal coupling. As shown in the bottom panel of Figure 2, the total energy is also almost conserved. The relative energy error is only 10⁻⁴. For the case with scattering opacity σ_s = 100, we find very similar behavior and accuracy. The only exception is that the gas temperature is unchanged and all the kinetic energy density is converted to the radiation energy density, as expected. Therefore, solutions from our transfer equation are consistent with the radiation moment equations. Notice that the last two terms in Equation (6), which we added in comparison to Mihalas & Klein (1982), are crucial to get the correct equilibrium state in a moving medium.

Zoom In Zoom Out Reset image size

Specific intensities not only allow calculation of the radiation energy density and flux, but also encode information about the angular distribution of the photons. An initially isotropic distribution of the specific intensity is shown in the left panel of Figure 3. After the equilibrium state is reached, the middle and right panels of Figure 3 show the spatial distributions of the specific intensities projected to the x − y plane for the cases with absorption opacity σ_a = 100 and scattering opacity σ_s = 100 respectively. The right-going intensities are increased while the left-going intensities are decreased. In other words, the specific intensities are beamed as expected.

Zoom In Zoom Out Reset image size

To check that the angular distribution of specific intensities agree with the Equations (17) and (25) quantitatively in the steady state, we set the left-hand side of the two equations to be zero and solve the two integral equations of the right-hand sides for all angles, which are shown as the black and blue ellipses in Figure 3. As this figure shows, the agreement is perfect.

The angular distribution of the photons can also be quantified by two components of the Eddington tensor f_xx ≡ P_{r, xx}/E_r and f_yy ≡ P_{r, yy}/E_r in this test. Figure 4 shows the evolution of f_xx and f_yy for the case with pure absorption opacity. They are 1/3 initially when the specific intensities are isotropic. In the equilibrium state, f_xx is increased to ∼0.37 while f_yy drops to ∼0.315, although the gas temperature is spatially isotropic. Consistent with Figure 3, the difference between f_xx and f_yy in the Eulerian frame is caused by the beaming of the fluid velocity. The degree of anisotropy is $\sim \mathcal {O}(v/{\mathbb {C}})^2$ in this thermal equilibrium case (Mihalas & Mihalas 1984).

Zoom In Zoom Out Reset image size

The steady state solution in a scattering opacity-dominated atmosphere is used to test the short characteristic module described in Davis et al. (2012), which solves the time-independent transfer equation and finds the solutions via iterations. This test is useful as it covers the transition from optically thick LTE to optically thin non-LTE regimes within the atmosphere. It also demonstrates that steady state solution can be achieved by solving the time-dependent RT equations.

The simulation domain is a cubic box with x × y × z ∈ (− 0.5, 0.5) × (− 0.5, 0.5) × (− 10, 10). The spatial resolution is fixed to be 16 × 16 × 1280 for the x, y, and z directions, respectively. Density of the atmosphere is chosen to be an exponential profile along vertical direction ρ(x, y, z) = 10⁻³exp (− z + 10). Gas temperature is fixed to be 1 everywhere and flow velocity is zero everywhere. Absorption opacity is σ_a = ρ while the scattering opacity is σ_s = (1 − )ρ. The parameter is the destruction parameter used in Davis et al. (2012), which is a constant here. If the Eddington tensor is fixed to be $1/3{\sf I}$, then the analytic solution for the radiation energy density is (Equation (30) of Davis et al. 2012)

$$\begin{eqnarray} &&E_r(\tau)=1-\frac{e ^{-\sqrt{3\epsilon }\tau }}{1+\sqrt{\epsilon }}, \end{eqnarray} \tag{ 33 }$$

where τ is the optical depth measured from the top z = 10. In order to be close to the Eddington approximation, we use one angle per octant. Periodic boundary conditions are used along the x and y directions. At the bottom of the simulation box, we copy the specific intensities along all angles from the last active zones to the ghost zones. At the top, only outgoing specific intensities are copied from the last active zones to the ghost zones while the incoming specific intensity is set to be zero in the ghost zones. The initial specific intensities are set to be 1/(4π) for all angles everywhere. The hydrodynamical equations are not evolved in this test. Because of the vacuum boundary condition at the top, the radiation energy density will decrease until cooling is balanced by the thermal emission in the equilibrium state. The smaller is, the smaller E_r is at the top. The solutions from our code at the equilibrium state for three different values of are shown as the solid black lines in Figure 5. Compared with the analytical solutions shown as the dashed red lines, they agree very well when changes from 0.1 to 10⁻⁴.

Zoom In Zoom Out Reset image size

Unlike the short characteristic module used in Davis et al. (2012), specific intensities are decomposed along each axis for the transport step. By propagating two beams of photons in vacuum, we demonstrate that this approach still allows accurate propagation of photons in a direction oblique to the grid.

We setup a 2D domain with box size L_x = 1, L_y = 4. Periodic boundary conditions are used along the x direction. Both absorption and scattering opacity are zero everywhere. We inject two beams from the bottom at 45 degrees with respect to the +x and −x axes, respectively. Because of the periodic boundary conditions, the two beams will cross each other and finally escape from the top. The distributions of E_r after the two beams have escaped for two different spatial resolutions are shown in Figure 6. Compared with Figure 6 of Davis et al. (2012), our approach gives similar results to the short characteristic method with quadratic interpolation. There is a small amount of numerical diffusion in the beam, although total flux remains conserved. The numerical diffusion is reduced with increasing resolution, which reflects convergence of the solution.

Zoom In Zoom Out Reset image size

As an aside, we note that M1 methods fail this test as the two beams of radiation merge into one at the point where they cross.

As described in the Appendix of Jiang et al. (2013b), the most useful test for our treatment of the transport step described in Section 4.2 is the dynamic diffusion test in an optically thick medium with pure scattering opacity. We setup a one-dimensional domain with L_x = 2. The density, temperature and flow velocity are all set to be 1. All the hydro quantities are uniform over the whole simulation domain and they do not evolve in the test. The dimensionless speed of light is ${\mathbb {C}}=10$. The value of ${\mathbb {P}}$ does not affect the test. The scattering opacity σ_s is a parameter that controls the diffusion coefficient $D\equiv {\mathbb {C}}/(3\sigma _s)$. We use 128 grid points and periodic boundary conditions.

The initial profiles of radiation energy density and flux for the region x ∈ (− 0.5, 0.5) are

$$\begin{eqnarray} E_r(x,0)&=&\exp\, (-40 x^2),\nonumber \\ F_r(x,0)&=&\frac{80 x}{3\sigma _s}E_r(x,0)+\frac{4v}{3{\mathbb {C}}} E_r(x,0). \end{eqnarray} \tag{ 34 }$$

For other positions, E_r and F_r are fixed to be exp (− 10) and $4vE_r/(3{\mathbb {C}})$. We will only focus on the evolution of the region $x\in (-0.5+\sqrt{Dt}, 0.5-\sqrt{Dt})$, as this part will not be affected by the finite domain of the simulation box. In order to be consistent with the Eddington approximation, so that analytical solutions can be used for comparison, we only use one angle per octant. Therefore the specific intensities are uniquely determined by E_r and F_r.

In the diffusion limit in the Eddington approximation, the profile of E_r at time t should be (Sekora & Stone 2010; Jiang et al. 2013b)

$$\begin{eqnarray} &&E_r(x,t)=\frac{1}{(160Dt+1)^{1/2}}\exp \left(\frac{-40(x-vt)^2}{160Dt+1}\right). \end{eqnarray} \tag{ 35 }$$

Comparisons between the numerical and analytical solutions at three different times for two different values of opacity σ_s are shown in Figure 7. This figure shows that our algorithm can calculate the diffusion and advection processes accurately for a wide range of opacity. In order to show the necessity of separating the diffusive and advective terms as described in Section 4.2, we repeat the test with optical depth per cell τ_s = 625 and use the original transport scheme described by Stone & Mihalas (1992). The result is shown in Figure 8. Compared with the left panel of Figure 7, the numerical diffusion rate is clearly dominant over the physical diffusion rate and our treatment shows significant improvement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to test the full radiation hydrodynamic algorithm, we carry out the linear wave test in a 2D domain based on the dispersion relation given by Jiang et al. (2012). We setup a uniform medium with background state ρ = T = P = 1 and γ = 5/3. Radiation energy density E_r = 1 and radiation flux is zero in the background state. In order to be consistent with the Eddington approximation used to calculate the dispersion relation, we use one angle per octant. The specific intensity in the background state is I = 1/(4π). We fix the dimensionless speed of light to be ${\mathbb {C}}=10$. The box size is L_x = 1 and L_y = 1 and the wave vector is aligned with x axis in order to calculate the propagation speed and damping rate more easily. As the hydrodynamics algorithm is second-order accurate, while the RT module is only first-order accurate due to operator splitting, we expect that the overall convergence rate for the whole code will depend on the parameter regime. For gas pressure dominated or optically thin regimes it will be second-order accurate, while for optically thick radiation pressure dominated regimes it will be first-order accurate. This is confirmed in Figure 9. We show the change of L1 error with spatial resolution for two sets of parameters $\tau _a=0.01, {\mathbb {P}}=1$ and $\tau _a=10, {\mathbb {P}}=10$. When the optical depth per wavelength is only 0.01, the code is second-order accurate. When the optical depth per wavelength is increased to 10 and radiation pressure is larger than gas pressure, the code is first-order accurate.

Zoom In Zoom Out Reset image size

In Figure 10, we compare the wave propagation speed and damping rate from our code with the analytical solutions as given in the Appendix B of Jiang et al. (2012) for optical depth per wavelength from 0.01 to 100 and ratio between radiation pressure and gas pressure from 0.01 to 100. In order to calculate the phase velocity and damping rate from our code, we evolve the linear wave for one period and calculate the Fourier transform of the density profile at the end of the simulation with background state subtracted. This is compared with the Fourier transform of the initial density profile with background state subtracted. We identify the positions of the component with maximum power in the phase space to be ϕ₀ and ϕ_t for the initial and final solutions. The values of maximum power are A₀ and A_t, respectively. Therefore, the difference between the phase velocity from the code and the analytical solution is δv = (ϕ_t − ϕ₀)/(kT), where k = 1/(2π) is the wave number and T is the period. The damping rate in the code can be calculated as ln (A₀/A_t)/T. The numerical values shown in Figure 10 are calculated with 512 grid points per wavelength. With this spatial resolution, our code can get the phase velocity and damping rate accurately over all the parameters we have explored. As ${\mathbb {C}}$ is only 10, which is the regime that the algorithm will be most useful, we can easily be in the regime where photons are well-coupled to the gas, and radiation pressure contributes a significant fraction of the restoring force for radiative acoustic modes. This is the case when τ_a > 10 with ${\mathbb {P}}=1$ and ${\mathbb {P}}=10$ in Figure 10, as the phase velocity is larger than the adiabatic sound speed for these cases. Again, our code captures these modes accurately.

Zoom In Zoom Out Reset image size

The effects of radiation on the structure of strong shocks have been used as standard tests for many radiation hydrodynamic codes (González et al. 2007; Sekora & Stone 2010; Zhang et al. 2011; Jiang et al. 2012). However, because of the complexity of the shock solutions, the numerical solutions are either compared with simple analytical estimates (Vaytet et al. 2013), or semi-analytical solutions based on the Eddington approximation in more quantitative tests (Lowrie & Edwards 2008; Jiang et al. 2012). Recently, McClarren & Drake (2010) have pointed out that the radiation field near the Zel'dovich spike in subcritical shocks can be very anisotropic, which invalidates the use of the Eddington approximation, and may provides a good diagnostic to test the accuracy of more advanced RT algorithms.

Recently, the steady-state structure of radiation modified shocks at a variety of Mach numbers has been worked out without any assumption regarding the Eddington tensor by Ferguson (2014). We adopt these new solutions to test our full RT and hydrodynamics algorithms.⁷

We initialize the shock solutions in one dimension with dimensionless pre-shock parameters ρ₀ = T₀ = E_{r, 0} = 1 and $v_0=\mathcal {M}$, where the Mach number $\mathcal {M}$ controls the relative importance of radiation to gas pressure. The other dimensionless numbers are ${\mathbb {C}}=1.73\times 10^3$, ${\mathbb {P}}=10^{-4}$, and σ_a = 577.4. These values are the same as in Section 6.2 of JSD12, so that the radiation shock solutions can be directly compared. Gas temperature T is calculated as T = P/(0.6ρ) and γ = 5/3. We use 4096 grid points and 10 angles per octant for all the solutions. All quantities on the left boundary are fixed to the pre-shock parameters, while on the right boundary all quantities in the ghost zones are copied from the last active zone. The steady state numerical solutions can be compared with the initial conditions to see how well the code can hold the semi-analytical solutions. Our numerical solutions after three flow crossing times for $\mathcal {M}=1.2,\ 2,\ {\rm and}\ 3$ are shown as black dots in Figures 11, 12, and 13 while the semi-analytical solutions are shown as dashed red lines in each figure. The numerical solutions agree with the semi-analytical solutions very well for the three cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the case with $\mathcal {M}=1.2$, gas temperature increases monotonically from upstream to downstream. The f_xx component of the Eddington tensor is close to 1/3 near the shock front, although it can be seen that f_xx can be smaller than 1/3 in the immediate post-shock flow, as found by Sincell et al. (1999) and JSD12. When the Mach number is increased to 2 and then 3, radiation pressure downstream becomes significantly stronger and the Zel'dovich spike (Zel'Dovich & Raizer 1967) emerges. In Figure 13 with $\mathcal {M}=3$, not only does the gas temperature T form a spike near the shock front, but also the radiation temperature T_r. At the same time, the radiation field becomes very anisotropic as the peak of f_xx differs from 1/3 by 50%. A similar spike in radiation temperature is evident in Figure 15 of JSD12 when the angular distributions of photons are calculated self-consistently with the VET approach. This feature is called "anti-diffusion" in McClarren & Drake (2010) because the direction of radiation flux at the shock front is along the same direction as the radiation energy density gradient. This is the opposite behavior to what is assumed in the diffusion approximation. It occurs because the RT equation requires ${{\boldsymbol F}}_r\propto -{{\boldsymbol \nabla}}\cdot {\sf P_r}=-{{\boldsymbol \nabla}}\cdot ({\sf f}E_r)$ in steady state. When the Eddington tensor has a strong spatial gradient as in the shock front, ${{\boldsymbol F}}_r$ is not necessarily in the same direction of $-{{\boldsymbol \nabla}}E_r$, as assumed in the diffusion approximation. These results point out that it is crucial to adopt algorithms that can treat the anisotropic nature of the radiation field accurately.

The anisotropic radiation field in our solution can also be demonstrated directly by plotting the specific intensities along different directions. Profiles of specific intensity along four different angles for the case with Mach number $\mathcal {M}=3$ are shown in Figure 14. The intensity along rays pointing in the −x direction is much larger than intensities in other directions in the upstream region. This is responsible for the pre-heating of the gas. The intensity along rays in the direction of the flow actually increase monotonically from upstream to downstream, which is quite different from the profile of T_r shown in Figure 13. The spike in radiation energy density near the shock front is dominated by the rays perpendicular to the flow direction. This is because just in front of the hot downstream gas, rays perpendicular to the flow direction have a longer path length through the narrow shell of hot gas, and therefore have proportionally more emitted photons than directions with shorter path lengths through the emission region. This is consistent with the fact that f_xx < 1/3 in Figure 13 near the shock front.

Zoom In Zoom Out Reset image size

The ability to capture shadows cast by an optically thick object in an optically thin environment requires propagating photons along different directions correctly. This test has been widely used to demonstrate the differences between FLD, M1, and RT algorithms based on short characteristics (Hayes & Norman 2003; González et al. 2007; McKinney et al. 2013, JSD12). Although the M1 method is able to capture the shadow formed by one beam, therefore demonstrating an advantage of M1 over FLD (González et al. 2007), it cannot propagate two beams correctly (McKinney et al. 2013). As our method directly solves for the specific intensity along different directions, it should be able to capture the shadow from multiple beams accurately, as was the case with the VET method in JSD12.

For this test, we set up a rectangular box of size (− 0.5, 0.5) × (− 0.3, 0.3) with dimensionless parameter ${\mathbb {C}}$ chosen to be 10 to speed up the calculation. The background medium has density ρ₀ = 1 and temperature T₀ = 1. Density inside the elliptical region r = x²/a² + y²/b² ⩽ 1 with a = 0.1 and b = 0.06 is ρ(r) = ρ₀ + (ρ₁ − ρ₀)/[1  +  exp  (10(r − 1))], where ρ₁ = 10ρ₀. Gas temperature is set to be 1/ρ. Absorption opacity σ_a = T^−3.5ρ² is used in this test. Beamed photons with intensity I = 103.13 are injected from the left x boundary at angles ±15° with respect to the x-axis. One beam at −15° and another beam at 15° with the same intensity are injected along the top and bottom y boundaries respectively. The injected beam angles are chosen to directly correspond to rays in our angular grid. We use n_a = 12 with all the angles uniformly distributed in the 2D plane as described at the end of Section 3. The intensities in the ghost zones along the right x boundaries are copied from the last active zones. We use a resolution of 512 × 256 cells for this test. The hydrodynamical equations are not evolved (although see Proga et al. (2014) for the hydrodynamical evolution of an irradiated cloud computed using our VET method). Although we use different values of ${\mathbb {C}}$ compared with the shadow test shown in Section 6.5 of JSD12, all the other important dimensionless numbers and angles of the incoming intensities are the same. For example, the photon mean free path inside the cloud is only 3.2 × 10⁻⁶ of the horizontal size of the simulation box due to the low temperature and high density, while it is the same as the box size outside the cloud. Therefore, results from the two different RT algorithms can be directly compared. Radiation energy density, xx and yy components of the Eddington tensor after the photons have propagated through the simulation box are shown in Figure 15. Compared with Figure 19 of JSD12, the umbra and penumbra calculated by the two codes are indeed very similar. Note also the results differ significantly from those shown by McKinney et al. (2013) computed using the M1 method. In order to see the effects of angular resolution, we have also done another test with the total number of angles increased to n_a = 36 while the angles of incoming beams remain the same. The result is almost identical. This is not surprising as the radiation field is dominated by the two incoming beams. However, the reradiated emission from the cloud, which is three orders of magnitude smaller than the incoming radiation field, is better resolved with more angles.

Zoom In Zoom Out Reset image size