A GODUNOV METHOD FOR MULTIDIMENSIONAL RADIATION MAGNETOHYDRODYNAMICS BASED ON A VARIABLE EDDINGTON TENSOR

To begin with, it is useful to summarize all of the individual steps in the overall algorithm.

Step 1. Using the gas variables to compute the opacities and source function, solve the transfer equation over a large set of angles using short characteristics, and compute the Eddington tensor ${\bm {\sf f}}$.

Step 2. Reconstruct the left and right states in the primitive variables at cell interfaces along each of the x-, y-, and z-directions independently. We reconstruct the primitive variables, instead of the conserved variables as in SS10, since we have found it to be less oscillatory and more accurate. Either second- or third-order spatial reconstruction is possible.

Step 3. Add the radiation energy and momentum source terms to the left and right states (see Equations (79) and (80) in SS10). Source terms for the left (right) states are calculated by using cell-centered quantities from the cell to the left (right) of the interface, and use the modified Godunov method.

Step 4. Compute 1D fluxes of the conserved variables with the appropriate Riemann solver. For radiation hydrodynamic simulations, we use the HLLC solver, while for radiation MHD simulations we generally use HLLD. We use the adiabatic sound speed (instead of the radiation-modified sound speed as was used by SS10) to calculate the fluxes. This adds extra dissipation and makes the algorithm more robust. We do not calculate fluxes for E_r and ${\mbox{\boldmath $F$}}_r$ because they are not updated in this step; therefore the Riemann solvers are the same as those described in Stone et al. (2008).

Step 5. For radiation MHD calculations, use CT algorithm, (Step 3 of the Athena 3D algorithm) to calculate the electric field at cell corners and update the face-centered magnetic field.

Step 6. Evolve the left and right states at each interface with the transverse flux gradients for a half time step Δt/2, as required for the unsplit CTU integrator in Athena.

Step 7. For radiation MHD calculations, calculate the cell-centered electric field at time t^{n + 1/2} to use as a reference state in CT algorithm (Step 6 of the Athena 3D algorithm).

Step 8. Calculate new fluxes along each direction with the corrected left and right states and the appropriate Riemann solver.

Step 9. Update the area-averaged magnetic field at cell faces from time step n to n + 1 using CT (Step 8 of Athena 3D algorithm).

Step 10. Update the gas variables from time n to n + 1 by adding the divergence of the flux gradient and the source terms using the modified Godunov method.

Step 11. Estimate the value of radiation energy source term added to the gas total energy using the method described in Section 3.3.1. Then update the radiation variables $E_r,\ {\mbox{\boldmath $F$}}_r$ using fully implicit, backward Euler differencing of the radiation moment Equation (8). This requires solving a large linear system in multidimensions. Gas variables do not change during this step.

Step 12. Update the time t^{n + 1} = tⁿ + Δt, calculate the new time step according to the CFL condition with the adiabatic sound speed (fast magnetosonic speed) for radiation hydrodynamics (MHD).

In comparison to the algorithm introduced in SS10, the primary changes we have made to extend the scheme to multidimensions are in Step 1 (compute the Eddington tensor in multidimensions), Step 2 (reconstruction in the primitive rather than conserved variables), Step 11 (using an estimate of the radiation energy source term, and solving the implicit moment equations in multidimension), and Step 12 (choose a time step based on the adiabatic sound speed).

In this subsection, we describe in detail the modified Godunov method to calculate the flux, in particular, how to add the stiff source terms to the left and right states.

3.2.1. Reconstruction Step

To compute the left and right states of the vector of conserved variables required to calculate the fluxes via a Riemann solver, we reconstruct the primitive variables as originally proposed by Miniati & Colella (2007), instead of the conserved variables (SS10). In hydrodynamics, the primitive and conserved variables are related via

$$\begin{eqnarray} &&\fontsize{7.5}{9.5}\selectfont\hbox{$\left(\begin{array}{@{}c@{}}\rho \\ v_x \\ v_y \\ v_z \\ P \end{array} \right)$}\nonumber\\ &&\quad \fontsize{7.5}{9.5}\selectfont\hbox{$ =\!\left[ \begin{array}{@{}c@{\hspace{7.3pt}}c@{\hspace{7.3pt}}c@{\hspace{7.3pt}}c@{\hspace{7.3pt}}c@{}}1 & 0 & 0 & 0 & 0\\ 0 & 1{/}\rho & 0 & 0 & 0 \\ 0 & 0 & 1{/}\rho & 0 & 0\\ 0 & 0 & 0 & 1{/}\rho & 0\\ (\gamma\,{-}\,1)\big(v_x^2+v_y^2+v_z^2\big)\!{\big/}\!2 & -(\gamma\,{-}\,1)v_x & -(\gamma\,{-}\,1)v_y & -(\gamma\,{-}\,1)v_z & \gamma\,{-}\,1 \end{array} \right]$}\nonumber\\ &&\qquad\fontsize{7.5}{9.5}\selectfont\hbox{$\times \left(\begin{array}{@{}c@{}}\rho \\ \rho v_x \\ \rho v_y \\ \rho v_z \\ E \end{array} \right).$} \end{eqnarray} \tag{ 10 }$$

The radiation source terms for velocity and pressure are

$$\begin{eqnarray} {\bf S}_{\rm r}({\mbox{\boldmath $v$}})&=&{-}{\mathbb {P}}{\bf S}_{\rm r}({\mbox{\boldmath $P$}})/\rho \\ \nonumber S_r(P) &=& (\gamma -1){\mathbb {P}}{\mbox{\boldmath $v$}}\cdot {\bf S}_{\rm r}({\mbox{\boldmath $P$}}) -(\gamma -1){\mathbb {P}}{\mathbb {C}}S_r(E), \end{eqnarray} \tag{ 11 }$$

where ${\bf S}_{\rm r}({\mbox{\boldmath $P$}})$ and S_r(E) are given by Equation (2). These source terms are the same in the case of radiation MHD, even though in this case the definition of total gas energy E includes the magnetic energy. Either second- or third-order reconstruction schemes can be used (Stone et al. 2008).

The radiation source terms for the primitive variables must be added to the left and right states for a half time step. This requires calculating the propagation operator $\newmathcal {I}$ to project off the unstable mode (e.g., Miniati & Colella 2007; SS10), using the gradient of radiation source terms on the plane of primitive variables ${\mbox{\boldmath $\nabla $}}_W{\bf S}_{\rm r}({\mbox{\boldmath $W$}})$. As discussed in SS10, in most cases it is the energy source term S_r(E) (or equivalently S_r(P)) that defines the stiffness of the problem. The momentum source term can be added as a normal body force. However, in the extreme case of radiation pressure so large that ${\mathbb {P}}\gtrsim {\mathbb {C}}$, the momentum source term may also become stiff. For example, this can happen in the inner regions of an accretion disk around a supermassive black hole, where ${\mathbb {C}}\sim 10^2$ while ${\mathbb {P}}\sim 10^2$ (e.g., Shakura & Sunyaev 1973; Turner et al. 2003). For this reason, we keep the leading term in the momentum source term, so ${\mbox{\boldmath $\nabla $}}_WS_r({\mbox{\boldmath $W$}})$ can be written as

$$\begin{equation} {\mbox{\boldmath $\nabla $}}_W{\bf S}_{\rm r}({\mbox{\boldmath $W$}})= \left(\begin{array}{@{}ccccc@{}}0 & 0 & 0 & 0 & 0\\ 0 & S^{v_x}_{v_x} & 0 & 0 & 0 \\ 0 & 0 & S^{v_y}_{v_y} & 0 & 0 \\ 0 & 0 & 0 & S^{v_z}_{v_z} & 0 \\ S^P_{\rho } & S^P_{v_x} & S^P_{v_y} & S^P_{v_z} & S^P_P \end{array} \right), \end{equation} \tag{ 12 }$$

where S^x_y ≡ ∂S_r(x)/∂y for any quantity x and y.

For example, in the x-direction $S^{P}_{v_x}/S^P_P\sim \newmathcal {O}(v/{\mathbb {C}}^2)$ and $S^{P}_{v_x}/S^P_{\rho }\sim \newmathcal {O}(v/{\mathbb {C}}^2)$. To order $v/{\mathbb {C}}$, we can take $S^P_{v_x}=S^P_{v_y}=S^P_{v_z}=0$, which significantly simplifies the analysis (a significant advantage of using the primitive variables). Then the propagation operator is

$$\begin{eqnarray} {\sf \newmathcal {I}}\left(\frac{\Delta t}{2}\right)&=&\frac{1}{\Delta t/2}\int _0^{\Delta t/2}e^{\tau {\mbox{\boldmath $\nabla $}}_W {\bf S}_{\rm r}({\mbox{\boldmath $W$}})}d\tau \nonumber \\ &=& \left(\begin{array}{@{}ccccc@{}}1 & 0 & 0 & 0 & 0 \\ 0 & \beta _x & 0 & 0 & 0 \\ 0 & 0 & \beta _y & 0 & 0 \\ 0 & 0 & 0 & \beta _z & 0 \\ (\alpha -1)S^P_{\rho }\big/S^P_P & 0 & 0 & 0 & \alpha \end{array} \right)\nonumber \\ &\approx &\left(\begin{array}{@{}ccccc@{}}1 & 0 & 0 & 0 & 0 \\ 0 & \beta _x & 0 & 0 & 0 \\ 0 & 0 & \beta _y & 0 & 0 \\ 0 & 0 & 0 & \beta _z & 0 \\ (1-\alpha)P/\rho & 0 & 0 & 0 & \alpha \end{array} \right), \end{eqnarray} \tag{ 13 }$$

where we define the following quantities:

$$\begin{eqnarray} \beta _x&=&\left[{\rm exp}\big(S^{v_x}_{v_x}\Delta t/2\big)-1\right]\big/\big(S^{v_x}_{v_x}\Delta t/2\big), \nonumber \\ \beta _y&=&\left[{\rm exp}\big(S^{v_y}_{v_y}\Delta t/2\big)-1\right]\big/\big(S^{v_y}_{v_y}\Delta t/2\big), \nonumber \\ \beta _z&=&\left[{\rm exp}\big(S^{v_z}_{v_z}\Delta t/2\big)-1\right]\big/\big(S^{v_z}_{v_z}\Delta t/2\big), \nonumber \\ \alpha &=&\left[{\rm exp}\big(S^P_P\Delta t/2\big)-1\right]\big/\big(S^P_P\Delta t/2\big). \end{eqnarray} \tag{ 14 }$$

Note that $S^{v_x}_{v_x}$, $S^{v_y}_{v_y}$, $S^{v_z}_{v_z}$$\sim \newmathcal {O}(-\sigma _t{\mathbb {P}}/{\mathbb {C}})$ and when $\sigma _t{\mathbb {P}}/{\mathbb {C}}\ll 1$, then β_x, β_y, β_z ≈ 1. In this case the propagation operator is reduced to the form used in SS10. When $\sigma _t{\mathbb {P}}/{\mathbb {C}}\gg 1$, then β_x, β_y, and β_z ensure that the algorithm is still stable in this regime. With this propagation operator, the source terms added to the left and right states are $0.5\Delta t\times {\sf \newmathcal {I}}(\Delta t/2){\bf S}_{\rm r}({\mbox{\boldmath $W$}})$.

3.2.2. Fluxes from the Riemann Solver

3.2.3. Updating the Conserved Variables with the Predict-correct Scheme

To achieve second-order accuracy, the conserved variables in this step are updated with a predict-correct scheme, similar to the approach taken by Miniati & Colella (2007) and SS10. The radiation quantities (energy and flux) have already been updated by the implicit solution of the radiation subsystem, so they are kept constant during this step. To be consistent with our reconstruction step, the stiffness of the momentum source terms in some regions must also be taken into account.

Given the divergence of the fluxes (computed as described above) ${\mbox{\boldmath $\nabla $}}\cdot {\mbox{\boldmath $F$}}^{n+1/2}$ and the radiation source term ${\bf S}_{\rm r}({\mbox{\boldmath $U$}}^{n})$, a predict solution is estimated as

$$\begin{eqnarray} &&\hat{{\mbox{\boldmath $U$}}}={\mbox{\boldmath $U$}}^n+\Delta t({\bm{\sf I}}-\Delta t {\mbox{\boldmath $\nabla $}}_U{\bf S}_{\rm r}({\mbox{\boldmath $U$}}^n))^{-1}({\bf S}_{\rm r}({\mbox{\boldmath $U$}}^n)-{\mbox{\boldmath $\nabla $}}\cdot {\mbox{\boldmath $F$}}^{n+1/2}),\nonumber\\ \end{eqnarray} \tag{ 15 }$$

where ${\mbox{\boldmath $\nabla $}}_U{\bf S}_{\rm r}({\mbox{\boldmath $U$}}^n)$ is the gradient of the radiation source term ${\bf S}_{\rm r}({\mbox{\boldmath $U$}})$ with respect to the conserved variables ${\mbox{\boldmath $U$}}^n$ as

$$\begin{eqnarray} {\mbox{\boldmath $\nabla $}}_U{\bf S}_{\rm r}({\mbox{\boldmath $U$}})= \left(\begin{array}{@{}ccccc@{}}0 & 0 & 0 & 0 & 0\\ 0 & S^{M_x}_{M_x} & 0 & 0 & 0 \\ 0 & 0 & S^{M_y}_{M_y} & 0 & 0 \\ 0 & 0 & 0 & S^{M_z}_{M_z} & 0 \\\ms S^E_{\rho } & S^E_{M_x} & S^E_{M_y} & S^E_{M_z} & S^E_E \end{array} \right). \end{eqnarray} \tag{ 16 }$$

The error in this predicted solution, (Δt), is estimated as

$$\begin{eqnarray} &&\epsilon (\Delta t)={\mbox{\boldmath $U$}}^n+\frac{\Delta t}{2}({\bf S}_{\rm r}({\mbox{\boldmath $U$}}^n) + {\bf S}_{\rm r}(\hat{{\mbox{\boldmath $U$}}})) -\Delta t {\mbox{\boldmath $\nabla $}}\cdot {\mbox{\boldmath $F$}}^{n+1/2}-\hat{{\mbox{\boldmath $U$}}}.\nonumber\\ \end{eqnarray} \tag{ 17 }$$

Then the correction step is

$$\begin{eqnarray} &&{\mbox{\boldmath $U$}}^{n+1}=\hat{{\mbox{\boldmath $U$}}}+({\bm{\sf I}}-\Delta t {\mbox{\boldmath $\nabla $}}_U{\bf S}_{\rm r}({\mbox{\boldmath $U$}}^n))^{-1}\epsilon (\Delta t). \end{eqnarray} \tag{ 18 }$$

At the end of this step, all gas and radiation quantities are updated to time step n + 1, and the entire algorithm can be repeated in the next cycle.

In this subsection, we describe the extensions to SS10 required to integrate the radiation subsystem Equation (8) in multidimensions.

3.3.1. Special Treatment of the Energy Source Term

The radiation energy density E_r and flux ${\mbox{\boldmath $F$}}_r$ are updated from time step n to n + 1 based on the updated gas quantities at time step n + 1, using first-order backward Euler time differencing to make the method stable with a time step much larger than the light crossing time. However, in SS10, the source terms on the right-hand side of Equation (8) were calculated using the temperature at time step n. We have found that this can introduce significant error in the total energy when the gas and radiation are far from thermal equilibrium, and the time step is much bigger than the thermalization time. The source of this error is the assumption that the gas temperature is constant during the update and thus the energy source terms added to the gas energy density and radiation energy density are not the same. This will be worse if there is continuous heating or cooling source terms and the energy error may accumulate. In reality, the gas temperature should evolve simultaneously with the radiation energy density until thermal equilibrium is reached.

To reduce this energy error, we first estimate the change of gas energy density due to the radiation energy source term in the modified Godunov step. Let ${\mbox{\boldmath $\nabla $}}\cdot {\mbox{\boldmath $F$}}^{n+1/2}$ be the Riemann flux calculated in Section 3.2.2. The updated gas total energy after the Godunov step is E^{n + 1} and the gas total energy at the beginning of the step is Eⁿ. Then the change of gas energy due to the radiation energy source term in this step dS_r(E) can be estimated as

$$\begin{eqnarray} &&dS_r(E)=E^{n+1}-(E^{n}-dt{\mbox{\boldmath $\nabla $}}\cdot {\mbox{\boldmath $F$}}^{n+1/2}). \end{eqnarray} \tag{ 19 }$$

To conserve total energy, the energy source term added to the radiation subsystem is then

$$\begin{eqnarray*} \fontsize{9.5}{11.5}\selectfont\hbox{$\tilde{S}_r(E)$}&\fontsize{9.5}{11.5}\selectfont\hbox{$\,{=}\,$}&\fontsize{9.5}{11.5}\selectfont\hbox{$\newmathcal {R}({\mathbb {C}}\sigma _a(\tilde{T}^4-E_r) +(\sigma _a-\sigma _s){\mbox{\boldmath $v$}}\cdot ({\mbox{\boldmath $F$}}_r-({\mbox{\boldmath $v$}}E_r+{\mbox{\boldmath $v$}}\cdot {\bm {\sf P}} _r)/\mathbb {C}))$}\nonumber\\ &&\fontsize{9.5}{11.5}\selectfont\hbox{$-\,(1-\newmathcal {R})dS_r(E)/(dt{\mathbb {P}}),$} \end{eqnarray*}$$

where the estimated temperature

$$\begin{eqnarray*} &&\tilde{T}^4=-(1/(dt{\mathbb {C}}\sigma _a)+1)dS_r(E)/{\mathbb {P}}+E_r^n. \end{eqnarray*}$$

Here $\newmathcal {R}$ is a parameter we can choose to get the best energy conservation. In practice, we find $\newmathcal {R}=0.05$ can give the best results for the tests we have done.

With this special treatment of the energy source term in the radiation moment equations, we can reduce the energy error significantly, especially for states that are initially far from thermal equilibrium in optically thick regions that are evolved with a time step much larger than the thermalization time (in practice, such states are extremely rare in any dynamical evolution, and usually occur only if set up specifically in the initial conditions). However, we still cannot conserve total energy to round-off error. We present tests in Section 4 that show in practice, the error in total energy conservation is small.

3.3.2. The 3D Matrix

The generalization of the fully implicit, backward Euler update of the radiation moment equations to 3D is straightforward. For this step, the vector of conserved quantities at time step n is Uⁿ_{i, j, k} = (E_r, F_{r, x}, F_{r, y}, F_{r, z}). The radiation flux F_{r, x}, F_{r, y}, F_{r, z} are advanced to time step n + 1 by solving

$$\begin{eqnarray} {\bf U}_{i,j,k}^{n+1}&=&{\bf U}_{i,j,k}^n-\frac{\Delta t}{\Delta x}\left[{\bf F}^{{\rm HLLE}}_{i+1/2,j,k}-{\bf F}^{{\rm HLLE}}_{i-1/2,j,k}\right] \nonumber\\ &&-\frac{\Delta t}{\Delta y}\left[{\bf F}^{{\rm HLLE}}_{i,j+1/2,k}-{\bf F}^{{\rm HLLE}}_{i,j-1/2,k}\right]\nonumber \\ &&-\frac{\Delta t}{\Delta z}\left[{\bf F}^{{\rm HLLE}}_{i,j,k+1/2}-{\bf F}^{{\rm HLLE}}_{i,j,k-1/2}\right] +\Delta t {\bf S}\big({\bf U}^{n+1}_{i,j,k}\big),\quad\ \ \end{eqnarray} \tag{ 20 }$$

while the radiation energy density E_r is updated as

$$\begin{eqnarray} E_{r,i,j,k}^{n+1}&=&E_{r,i,j,k}^n-\frac{\Delta t}{\Delta x}\left[{\bf F}^{{\rm HLLE}}_{i+1/2,j,k}-{\bf F}^{{\rm HLLE}}_{i-1/2,j,k}\right]\nonumber\\ && -\frac{\Delta t}{\Delta y}\left[{\bf F}^{{\rm HLLE}}_{i,j+1/2,k}-{\bf F}^{{\rm HLLE}}_{i,j-1/2,k}\right]\nonumber \\ &&-\frac{\Delta t}{\Delta z}\left[{\bf F}^{{\rm HLLE}}_{i,j,k+1/2}-{\bf F}^{{\rm HLLE}}_{i,j,k-1/2}\right] +\tilde{S}_r(E), \end{eqnarray} \tag{ 21 }$$

where the F^HLLE are the vector of fluxes for the conserved quantities at each cell interface computed by an HLLE Riemann solver (Equation (39) in SS10) and dS_r(E) is the estimate energy source term (Equation (19)). Since the backward Euler differencing is only first order in time, first-order spatial reconstruction is used to compute the left and right states for the HLLE solver. Moreover, since both the HLLE fluxes and the source term S(U^{n + 1}) in Equations (20) and (21) are linear in the unknowns U^{n + 1}_{i, j, k}, then only one matrix solve is required per time step, which can represent a significant saving compared to the implicit solution of nonlinear equations that can result from other splittings (e.g., Turner & Stone 2001; González et al. 2007). The matrix of coefficients that must be inverted to solve the linear system in this step is given in Appendix A.

In order to calculate the VET, we use a formal solution of time-independent transfer Equation (6). The methods for solving the RT equation are described in Davis et al. (2012), and the reader is referred to that work for further details. At each time step, and for each grid cell, the specific intensity I_r must be integrated along many different rays $\hat{n}$. The angular discretization and corresponding quadratures are chosen to cover the unit sphere as uniformly as possible, but still be invariant under 90° rotations about the coordinate axes (e.g., Bruls et al. 1999). Along each ray, the method of short characteristics (e.g., Mihalas et al. 1978; Olson & Kunasz 1987) is used to integrate the transfer equation across the entire mesh. For multidimensional problems, this requires the interpolation of intensities, opacities, and emissivities from (only) neighboring grid zones. Simulations with scattering opacity (non-LTE problems) are solved iteratively with accelerated lambda iteration, with methods similar to those described in Trujillo Bueno & Fabiani Bendicho (1995). Intensity boundary conditions and parallelization are handled via ghost zones, as described in Section 3.5 of Davis et al. (2012).

During the integration along each angle, the zeroth moment (J) and second moment (${\bm {\sf K}}$) of the specific intensity are summed into running quadratures. This saves having to allocate an array to store I_r over all angles. For each discreet ray $\hat{n}_k$, there is a vector of direction cosines (μ_0k, μ_1k, μ_2k) with $\mu _{ik}=\hat{n}_k \cdot \hat{x}_i$ and quadrature weights w_k. The moments are then given by

$$\begin{eqnarray} J &=& \sum _{k=0}^{n_a-1} w_k I_k\qquad \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} K_{ij} &=& \sum _{k=0}^{n_a-1} w_k I_k \mu _{ik} \mu _{jk}, \end{eqnarray} \tag{ 23 }$$

where $I_k \equiv I_{r}(\hat{n}_k)$. Up to a common factor of 4π/c, J and ${\bm {\sf K}}$ are equivalent to radiation energy density and pressure, respectively. Since they are computed by the radiation transfer solver, they will (in general) differ from the E_r and ${\bm {\sf P}}_r$ of the integrator. As defined in Equation (6), the VET is then simply the ratio of these quadratures ${\bm {\sf f}} = {\bm {\sf K}}/J$, evaluated in each grid cell.

Evolution of a Marshak wave is a 1D non-equilibrium diffusion process originally proposed by Marshak (1958). A semi-analytic solution is given by Su & Olson (1996), and we compare the results from our code with this solution. This test has also been used by many other authors (González et al. 2007; Gittings et al. 2008; van der Holst et al. 2011; Zhang et al. 2011). It consists of a cold uniform medium with constant absorption opacity σ_a. A constant radiation flux F^inc_r is applied at the left boundary at t = 0 and allowed to diffuse through the medium. The cold gas is heated up by the radiation field, and eventually the gas and radiation field reach thermal equilibrium. This is not a dynamical test, so the hydrodynamic algorithm is not used, and only the radiation subsystem is solved, with the gas temperature updated according to Equation (124) of SS10. All of the other parameters are the same as in SS10. In particular, the time step is limited by the light crossing time of each cell in a domain of size L = 0.1 and 512 grid points. As shown in Figure 3, our numerical solution matches the semi-analytic solution extremely well.

Zoom In Zoom Out Reset image size

The "tophat" or "crooked pipe" test (e.g., Gentile 2001) is designed to study the propagation of a radiation field in a low-opacity pipe with a complex shape. Surrounding the pipe is a high-opacity material. It is a challenging problem because it requires following the propagation of radiation through complex shapes bounded by discontinuities in opacities. We find this test useful to show that our first-order backward Euler differencing with VET can capture sharp gradients in the radiation field, as well as follow the propagation of radiation along the pipe properly.

The problem is initialized following the description in Gentile (2001), except that we solve the problem in Cartesian instead of cylindrical coordinates. Because of this difference, we cannot reproduce the solution given by Gentile (2001) exactly. Nonetheless, the test is still useful to demonstrate that we achieve qualitatively the same solution, and that the VET method can solve this problem properly.

The size of the simulation domain is (0, 7) × (− 2, 2) with a resolution of 512 × 256 grid points. Dense, opaque material with density 10 g cm⁻³ and opacity σ_a = 2000 cm⁻¹ is located in the following regions: (3, 4) × (− 1, 1), (0, 2.5) × (− 2, −0.5), (0, 2.5) × (0.5, 2), (4.5, 7) × (− 2, −0.5), (4.5, 7) × (0.5, 2), (2.5, 4.5) × (− 2, −1.5), and (2.5, 4.5) × (1.5, 2.5). The pipe, with density 0.01 g cm⁻³ and opacity σ_a = 0.2 cm⁻¹, occupies all other regions. The structure of the pipe is shown by the black line in the top panel of Figure 4. The dimensionless speed of light is 300. Initially, the material has a temperature 0.05 keV everywhere, and the radiation and material temperature are in equilibrium. A heating source with a fixed temperature 0.5 keV is located on the left boundary for −0.5 < y < 0.5. All other boundary conditions are outflow. We use the short-characteristic module to calculate the VET. An isotropic incoming specific intensity is applied only along the left boundary in the region covering the heating source. The incoming specific intensity is zero for all other boundaries. During the evolution, velocities are always kept zero so the material density does not change. Thus, the modified Godunov step to update the material quantities is not needed in this test, and we simply evolve the material temperature through the following two equations

$$\begin{eqnarray} c_v\rho (T^{n+1}-T^{n})&=&{-}dt{\mathbb {C}}\sigma _a\big(a_rT^{n+1,4}-E_r^{n+1}\big), \nonumber \\ E_r^{n+1}-E_r^{n}&=&dt{\mathbb {C}}\sigma _a\big(a_rT^{n+1,4}-E_r^{n+1}\big), \end{eqnarray} \tag{ 24 }$$

where the heat capacity c_v = 10¹⁵ erg g⁻¹ keV⁻¹. Note that we only update material temperature in this step; the radiation energy density is unchanged. The radiation energy source term added in this step is dS_r(E) = c_vρ(T^{n + 1} − Tⁿ). The radiation energy density and flux are then evolved using our first-order backward Euler update. We call this algorithm method I.

Zoom In Zoom Out Reset image size

First, we fix the time step to be 4 × 10⁻⁵ to resolve the light crossing time. Snapshots of the temperature distribution at two different times are shown in Figure 4. Note that the shadows around the corners of the pipe are captured properly. Five probes are placed at (x = 0.25, y = 0), (x = 2.75, y = 0), (x = 3.5, y = 1.25), (x = 4.25, y = 0), and (x = 6.75, y = 0) to monitor the change of the temperature in the pipe. The time evolution of the temperatures at the five points is shown as solid black lines in Figure 5. Compared with Figure 9 of Gentile (2001), which shows the result calculated with a Monte Carlo method, our first-order backward Euler method with VET gives very similar results when the time step is small enough to resolve the light crossing time. For the fifth point, note that it cools off slightly before being heated by the radiation wave, in agreement with the behavior noted by Gentile (2001). The most obvious difference in our solution compared to that of Gentile (2001) is that the temperature of the first two points increases too quickly. This is likely due to the fact that we neglect the light travel time in our short-characteristic module when we calculate the VET. Thus the first two points are affected by the heating source too early. We also find the evolution of the last three points to differ in detail from that shown by Gentile (2001); however for these points the difference between propagation in cylindrical versus planar (Cartesian) geometry may be important. In our solution, there is no geometrical dilution of the flux as it propagates radially, and therefore the heating rate per unit surface area on the walls of the pipe is increased. This can affect the detailed temporal evolution at the third through fifth probes. Most importantly, note that from Figure 4 our first-order implicit update can maintain very sharp gradients in the temperature at the walls of the pipe.

Zoom In Zoom Out Reset image size

As a second test of our method, we repeat the problem using a time step of 10⁻³ initially, and increasing it by a factor of 1.1 at each step until it reaches 1. The result is shown as the dashed line in Figure 5. With this larger time step, the error in the temporal evolution of the probes is increased, although the basic evolution is still correct. This is consistent with our expectations. If we are not interested following evolution on a light crossing or thermalization time, we can use a large time step and still maintain a stable solution with the correct asymptotic behavior. Instead, if we want to explore the evolution on these timescales, we must use a small enough time step to resolve them.

A second method to evolve the material temperature with our algorithms is to add the radiation energy source term based on the solution of the time-independent transfer equation used to calculate the VET. Details of how to compute the material heating rate from the solution of the transfer equation are given in Davis et al. (2012). We call this algorithm method II. The solution computed using this method for the tophat test is shown as red lines in Figure 5. The time step is fixed to be 10⁻³. The result does not change substantially when we decrease the time step further. This method also obtains a very similar solution to that shown in Figure 9 of Gentile (2001). Once again, the material is heated up too quickly everywhere but the location of the first probe. Again, this is likely due to the fact that the short-characteristic module solves the time-independent radiation transfer equation, so that we ignore the delay due to the finite speed of propagation of radiation down the pipe. The last three points also differ because of our use of a planar geometry. Nevertheless, method II still captures the basic evolution correctly.

Other stringent test of our radiation transport module used to calculate the VET (such as the beam test), are described in detail in Davis et al. (2012), and will not be repeated here.

The dispersion relation for linear wave solutions to the radiation hydrodynamic equations in a uniform homogeneous background medium, and assuming the Eddington approximation (fixed Eddington factor ${\sf f}=1/3{\sf I}$), have been analyzed by a variety of authors (e.g., Mihalas & Mihalas 1984; Bogdan et al. 1996; Johnson & Klein 2010). These authors solve the boundary value problem in which a complex wavenumber k is found for each real frequency ω. Johnson & Klein (2010) have used these solutions to test the Lagrangian radiation hydrodynamic code KULL by driving the waves with a time-dependent boundary condition; similar tests were used earlier by Turner & Stone (2001) to test the FLD module in the ZEUS code.

In contrast, Lowrie et al. (1999) give the dispersion relation for linear waves (again in the Eddington approximation) as an initial value problem, that is solving for the imaginary frequency ω at fixed real k. Although their analysis focused on the properties of linear waves in a moving fluid (they point out that in radiation hydrodynamics, linear wave solutions are not Galilean invariant), their approach is ideally suited for code tests. It allows the propagation of linear waves to be followed in a periodic domain, free from the complexity of the implementation of driving boundary conditions in an Eulerian code, which may in and of itself introduce spurious error into the solutions. The dispersion relation for linear waves in radiation MHD has also been considered by Blaes & Socrates (2001, 2003). We use both the hydrodynamic and MHD linear wave solutions for convergence testing in one, two, and three dimensions in the following subsections. We emphasize that in order to keep the solutions to the dispersion relation analytic so that they can be used in the test problems, the Eddington approximation must be adopted. There may be regions of parameter space where the Eddington approximation is not appropriate, and the properties of linear waves in this regime should be computed using a variable Eddington factor computed self-consistently with flow. However, such solutions are beyond the scope of this paper.

It is not possible to simply initialize this test problem using previous solutions for linear waves given in the literature, because there may be differences in the frame of reference, and in the order of the source terms kept, in the fundamental system of dynamical equations (our Equation (1)). For this reason, we have rederived the dispersion relation for linear waves in the system of equations solved by our code. Even with the Eddington approximation, these solutions are non-trivial, requiring finding the roots of a fourth-order polynomial for complex k in hydrodynamics, and a tenth-order polynomial for k in MHD. Appendices B and C give the details of our solutions in the case of hydrodynamics and MHD, respectively, while in Tables 1–3 we give the complete eigenvector for linear waves for several parameter values for use by others.

Table 1. Examples of Solutions to the Dispersion Relation (B1)

${\mathbb {P}}$	0.01	0.01
σ	0.01	10
δρ	10−3	10−3
δv	$\hphantom{-}1.27177\times 10^{-3}+8.15409\times 10^{-5} i$	1.00012 × 10−3 + 6.91295 × 10−6i
δP	$\hphantom{-}1.61075\times 10^{-3}+2.07402\times 10^{-4} i$	1.00019 × 10−3 + 1.36663 × 10−5i
δEr	$\hphantom{-}1.79137\times 10^{-8}+8.56498\times 10^{-9} i$	6.64619 × 10−7 + 4.83820 × 10−5i
δFr	−1.32035 × 10−6 + 3.88814 × 10−6i	−9.99976 × 10−6 + 1.40748 × 10−7i
ω	7.99077 + 0.512336i	6.28393 + 4.34354 × 10−2i
${\mathbb {P}}$	100	100
σ	0.01	10
δρ	10−3	10−3
δv	$\hphantom{-}1.00000\times 10^{-3}+8.93099\times 10^{-8} i$	9.99947 × 10−4 + 1.07703 × 10−5i
δP	$\hphantom{-}1.00000\times 10^{-3}+1.57240\times 10^{-7} i$	9.99998 × 10−4 + 1.60499 × 10−7i
δEr	−3.55851 × 10−13 + 6.41394 × 10−10i	−6.60096 × 10−9 + 6.41367 × 10−7i
δFr	−1.00000 × 10−9 + 2.10690 × 10−13i	−1.00130 × 10−9 + 5.35966 × 10−11i
ω	6.28319 + 5.61151 × 10−4i	6.28285 + 6.76716 × 10−2i
${\mathbb {P}}$	100	1
σ	0.1	10
δρ	10−3	10−3
δv	$\hphantom{-}1.00000\times 10^{-3}+1.15555\times 10^{-7} i$	1.00000 × 10−3 + 3.32795 × 10−7i
δP	$\hphantom{-}1.00000\times 10^{-3}+1.73114\times 10^{-8} i$	1.00000 × 10−3 + 2.94229 × 10−7i
δEr	−1.0093 × 10−12 + 6.41394 × 10−9i	3.419000 × 10−10 + 1.11408 × 10−6i
δFr	−1.00000 × 10−9 + 5.51807 × 10−13i	−1.00000 × 10−7 + 1.22263 × 10−10i
ω	6.28319 + 7.26052 × 10−4i	6.28319 + 2.09101 × 10−3i

Download table as:  ASCIITypeset image

6.1.1. Linear Waves in Radiation Hydrodynamics

We begin by studying the propagation of linear waves in 1D with hydrodynamics. We consider a uniform, homogeneous medium with background state ρ = T = P = 1, adiabatic index γ = 5/3, and zero velocity initially. The radiation flux is zero, and the radiation energy density E_r = 1 initially for all of the calculations. We assume a constant absorption opacity σ_a (perturbations of the opacity will always be second order, which can be neglected), and zero scattering opacity. In order to adopt the Eddington approximation, we fix the Eddington tensor to be diagonal with components f_ii = 1/3. The dimensionless speed of light ${\mathbb {C}}=10^4$. We use a domain of size L = 1 with 512 grid zones, and periodic boundary conditions for all variables at each edge. We initialize a wave solution (Appendix B) by adding an eigenmode (see Table 1) with a wavelength equal to the size of the domain L, and an amplitude of 10⁻⁶.

We perform a series of calculations in which we study the effect of varying the radiation to gas pressure ratio by varying the dimensionless pressure ${\mathbb {P}}$, while keeping the initial conditions fixed as above. We also study the effect of varying the optical depth per wavelength τ_a = σ_aL by varying the absorption opacity σ_a. For each of these calculations, we measure the phase velocity of the waves by determining how long it takes a fixed feature in the wave (the density maximum) to propagate a distance equal to one wavelength (L = 1). We also measure the damping rate of the wave by determining the best-fit amplitude to the entire wave profile after one period. We then compare these results to linear theory.

The properties of radiation-modified linear acoustic waves vary dramatically with optical depth and pressure ratio. Two dimensionless numbers can be used to define various regimes. The first, ${\mathbb {P}}{\mathbb {C}}\tau _a$, measure the importance of energy exchange between the radiation field and the material, while the second, ${\mathbb {P}}\tau _a$, measures the importance of momentum exchange. When both ${\mathbb {P}}{\mathbb {C}}\tau _a \ll 1$ and ${\mathbb {P}}\tau _a \ll 1$, energy exchange between the radiation field and material is very weak, the momentum of the radiation field can be neglected, so the waves are weakly damped and propagate at nearly the adiabatic sound speed. The damping rate increases with optical depth until τ_a ∼ 1; thereafter the radiation and material energy densities are strongly coupled, and the damping rate decreases with increasing optical depth. When ${\mathbb {P}}{\mathbb {C}}\tau _a \gg 1$ and ${\mathbb {P}}< 1$, the radiation and material energy densities are strongly coupled but the radiation momentum is still negligible, therefore in this case the damping rate decreases with optical depth until τ_a ∼ 1, and increases thereafter. Finally, if ${\mathbb {P}}{\mathbb {C}}\tau _a \gg 1$ and ${\mathbb {P}}> 1$, the momentum in the photons becomes important, and the damping rate reaches a minimum when ${\mathbb {P}}\tau _a\sim 1$ and increase afterward. To ensure our algorithms are accurate over a wide range of regimes, it is important to perform tests that span these dimensionless parameter values. These dimensionless quantities also clarify a potential shortcoming with numerical algorithm based on the reduced speed of light approximation (e.g., Gnedin & Abel 2001). Any arbitrary reduction in the speed of light will reduce ${\mathbb {C}}$ and the above dimensionless numbers correspondingly, which will alter the damping rate and phase velocity of the radiation-modified acoustic wave.

In addition to ${\mathbb {P}}$ and ${\mathbb {C}}$, the Boltzmann number Bo is another dimensionless number which is useful to quantify the relative importance of radiative and material energy transport in a radiating flow (e.g., Mihalas & Mihalas 1984); it is defined to be the ratio between the material enthalpy flux and radiative flux from a free surface. It is related to our dimensionless numbers ${\mathbb {P}}$ and ${\mathbb {C}}$ as

$$\begin{eqnarray} &&{\rm Bo}=\frac{4\tilde{v}\gamma }{\mathbb {PC}(\gamma -1)}, \end{eqnarray} \tag{ 25 }$$

where $\tilde{v}\equiv v/a_0$ is the typical flow velocity in units of a₀. In the linear wave tests, the typical flow velocity v can be replaced with the adiabatic sound speed. When Bo is small, energy transport is dominated by radiation field and material energy transport is dominant when Bo is large.

Figure 6 compares the numerically measured phase velocity and damping rate (shown as stars) for two different radiation pressures, and over a wide range of optical depths τ_a, in comparison to the solution of the linear dispersion relation (Equation (B1)) shown as solid lines. The parameters are chosen such that for the runs with ${\mathbb {P}}= 10^{-4}$, the dimensionless number ${\mathbb {P}}{\mathbb {C}}\tau _a$ spans $10^{-2} \le {\mathbb {P}}{\mathbb {C}}\tau _a \le 10^{2}$, while for the runs with ${\mathbb {P}}= 1$ both of these limits are a factor of 10⁴ larger. Based on the discussion above, for ${\mathbb {P}}= 10^{-4}$ we expect the damping rate to increase until τ_a ∼ 1 and decrease thereafter, while for ${\mathbb {P}}= 1$ we expect the opposite. The solid lines from the analytic solution of the dispersion relation clearly show these trends. In addition, the numerically measured phase velocity agrees with the analytic results over the entire range of optical depths for both values of the radiation pressure, while the damping rate is accurate in all cases except when ${\mathbb {P}}=1$ and τ_a ⩾ 1. However, for these parameter values, the damping rate is small and the numerical measured damping rates are not converged at the resolution used for this plot. We have confirmed that if we increase the resolution, the numerical measured damping rate converges to the analytic values. At the same time, we have also found when τ_a ⩾ 100, our algorithm no longer reproduces the correct damping rate. We discuss this further, and suggest a solution, below. When ${\mathbb {P}}=10^{-4}$, Bo = 10 and it is 10⁻³ when ${\mathbb {P}}=1$. In the case ${\mathbb {P}}=100$, Bo is only to 10⁻⁵.

Zoom In Zoom Out Reset image size

In Figure 7 we test linear wave propagation over a wide range of radiation pressures. Once again, the left panel compares the numerically measured phase velocity and damping rate (shown as stars) for three different radiation pressures, and over a wide range of optical depths τ_a, in comparison to the solution of the linear dispersion relation (Equation (B1)) shown as solid lines. Even though the properties of radiation-modified linear acoustic waves vary dramatically over this range of optical depths and pressure ratios, there is good agreement between the numerical and analytic solutions, except when ${\mathbb {P}}=1$ and τ_a > 1. For these parameters, we also find at higher resolution, our algorithm converges to the analytic damping rate, suggesting that the discrepancy is dominated by numerical diffusion at this resolution.

Zoom In Zoom Out Reset image size

Note that in Figure 7 we do not plot the numerical solution for τ = 100 and ${\mathbb {P}}=100$ or 0.01. We have found that the damping rate of linear waves stops converging for these parameter values. Our tests reveal that the problem is associated with our implicit solution of the radiation subsystem at a time step determined by the adiabatic sound speed in the material. In this regime, an accurate solution requires a time step which resolves the light crossing time across a cell. Note we find that this is only an issue for very large absorption opacity: large scattering opacities do not couple the radiation and material energy densities. When the optical depth per wavelength due to the absorption opacity becomes so large (${\mathbb {P}}\tau _a\sim {\mathbb {C}}$) that the radiation field and the material can be treated like a single fluid, we have found that it is more accurate to solve a single system of conservations laws for the total (radiation plus material) energy and momentum, rather than splitting the radiation subsystem from the material conservation laws.

As our next test, we measure the convergence rate of linear waves in 3D, for two different values of the problem parameters. We use a computational domain of size 2L × L × L, with the same number of grid cells per unit length L in each direction. We initialize the 1D solution inclined to the grid using the same technique to compute the appropriate volume averages of the solution as described in Stone et al. (2008). We propagate the wave for one period, and then measure the L1 error in the solution from

$$\begin{equation} \delta {\bf q} = \frac{1}{N} \sum _{i} \big\vert {\bf q}_{i} - {\bf q}_{i}^{0} \big\vert, \end{equation} \tag{ 26 }$$

where q⁰_i is the initial solution, and the sum is taken over all grid points. We repeat this calculation for a variety of different numerical resolutions (grid cells per unit length L), and plot the change in the L1 error with resolution. The result is shown in the right panel of Figure 7 for a radiation-dominated, optically thick fluid (${\mathbb {P}}=10^{2}$, τ_a = 10), denoted by the solid circles, and for a gas-pressure-dominated, optically thin fluid (${\mathbb {P}}=10^{-2}$, τ_a = 10⁻²), denoted by the open squares. In the radiation-dominated case the errors converge close to first order, while in the gas-pressure-dominated case, the errors converge close to second order. This behavior is expected. The errors in the former are dominated by the solution to the radiation subsystem, which uses a first-order accurate backward Euler step for stability. Thus, we expect the overall rate of convergence should be first order. On the other hand, in the gas-pressure-dominated case, the errors are dominated by the solution to the material conservation laws, which uses a second-order accurate modified Godunov step. In this case, the overall rate of convergence should be second order. Of course, it might be possible to increase the rate of convergence in the implicit solution of the radiation moment equations by adopting higher-order implicit differencing. However, limiting such methods to enforce monotonicity can be problematic (e.g., SS10). We prefer to adopt unconditionally stable first-order methods for the implicit solver.

Linear wave solutions when only radiation energy source terms are added with radiation momentum source terms neglected are also given by Davis et al. (2012) and used to test our radiation transfer module (see Figure 9 of that paper). Because the radiation energy source term is added explicitly, the radiation diffusion mode needs to be resolved to make the code stable, which can limit the time step significantly when Bo ≲ 1. Therefore the algorithm of Davis et al. (2012) is most suitable for the regime Bo > 0.01. Because we include all the radiation momentum and energy source terms, the dispersion relations given in this paper differ from Figure 8 of Davis et al. (2012). However, in regimes of parameter spaces which overlap, we do get very similar results, for example when ${\mathbb {P}}=10^{-4}$ and Bo = 10 shown in Figure 6.

6.1.2. Linear Waves in Radiation MHD

Next, in order to test the accuracy of our MHD algorithms with radiation, we consider the propagation of linear modes of radiation-modified magnetosonic waves. We do not consider Alfvén waves in this subsection since they are incompressible and are affected less by radiation than magnetosonic modes. In the case of radiation-modified MHD waves, the transverse components of the velocity must be included even for 1D solutions. However, as implemented our code does not include the transverse velocities in 1D during the implicit solution of the radiation moment equations. Thus, all of the MHD tests presented here have been performed in 2D using a grid of either 512 × 64, or in some cases 1024 × 64 cells.

For these MHD tests, we use a uniform, homogeneous medium with the hydrodynamic and radiation variables set to the identical values as were used for the hydrodynamic test (see the previous subsection). The strength of the magnetic field is $B_{0} = \sqrt{5/3}$ (which gives a ratio of the Alfvén to sound speed of 2). The direction of the magnetic field is 45° to the x-axis, and it is confined to the x–y plane (so B_z = 0). We use a 2D domain of size L × L with L = 1 with periodic boundary conditions for all variables at each edge. As before, we initialize a wave solution (Appendix C) by adding an eigenmode (see Tables 2 or 3) propagating in the x-direction with a wavelength equal to then size of the domain L, and an amplitude of 10⁻⁶. Once again, we perform a series of calculations in which we study the effect of varying the radiation to gas pressure ratio by varying the dimensionless pressure ${\mathbb {P}}$, and the effect of varying the optical depth per wavelength τ_a = σ_aL by varying the absorption opacity σ_a.

Table 2. Examples of Solutions for Radiation-modified Slow Mode

${\mathbb {P}}$	0.01	0.01
σ	0.01	10
δρ	10−3	10−3
δvx	7.77601 × 10−4 + 5.56463 × 10−5i	6.53191 × 10−4 + 2.46133 × 10−6i
δvy	1.20178 × 10−3 + 1.84722 × 10−4i	8.77920 × 10−4 + 5.58924 × 10−6i
δP	1.52579 × 10−3 + 2.97056 × 10−4i	1.00009 × 10−3 + 8.92617 × 10−6i
δBy	−7.15903 × 10−4 − 1.63064 × 10−4i	−4.44181 × 10−4 − 4.50308 × 10−6i
δEr	1.54496 × 10−8 + 1.02737 × 10−8i	3.069869 × 10−7 + 3.16009 × 10−5i
δFr, x	−1.89110 × 10−6 + 3.34728 × 10−6i	−6.53138 × 10−6 + 6.48914 × 10−8i
δFr, y	1.61623 × 10−7 + 1.67917 × 10−8i	1.17056 × 10−7 + 7.40429 × 10−10i
ω	4.88581 + 0.349636i	4.10401 + 1.54776 × 10−2i
${\mathbb {P}}$	100	100
σ	0.01	10
δρ	10−3	10−3
δvx	$\hphantom{-}6.53158\times 10^{-4}+3.18349\times 10^{-8} i$	$\hphantom{-}6.53158\times 10^{-4}+3.83780\times 10^{-6}i$
δvy	$\hphantom{-}8.77864\times 10^{-4}+7.22107\times 10^{-8} i$	$\hphantom{-}8.77819\times 10^{-4}+8.70693\times 10^{-6} i$
δP	$\hphantom{-}1.00000\times 10^{-3}+1.02703\times 10^{-7} i$	$\hphantom{-}9.99999\times 10^{-4}+1.04836\times 10^{-7}i$
δBy	−4.44142 × 10−4 − 5.81571 × 10−8i	−4.44093  ×  10−4 − 7.01467  ×  10−6i
δEr	−1.47792 × 10−13 + 4.18932 × 10−10i	−2.33054 × 10−9 + 4.18933 × 10−7i
δFr, x	−6.53158 × 10−10 + 9.61494 × 10−14i	−6.53472  ×  10−10 + 2.36270  ×  10−11i
δFr, y	$\hphantom{-}1.16852\times 10^{-7}-4.78590\times 10^{-9} i$	$\hphantom{.}1.17043\times 10^{-7}+1.15612\times 10^{-9}i$
ω	4.10391 + 2.000251 × 10−4i	4.10391 + 2.41136 × 10−2i
${\mathbb {P}}$	100	1
σ	0.1	10
δρ	10−3	10−3
δvx	$\hphantom{-}6.53158\times 10^{-4}+4.11763\times 10^{-8} i$	$\hphantom{-}6.53158\times 10^{-4}+1.18584\times 10^{-7}i$
δvy	$\hphantom{-}8.77864\times 10^{-4}+9.34194\times 10^{-8} i$	$\hphantom{-}8.77865\times 10^{-4}+2.69043\times 10^{-7} i$
δP	$\hphantom{-}1.00000\times 10^{-3}+1.13071\times 10^{-8} i$	$\hphantom{-}1.00000\times 10^{-3}+1.92178\times 10^{-7}i$
δBy	−4.44142  ×  10−4 − 7.52612  ×  10−8i	−4.44142 × 10−4 − 2.16754 × 10−7i
δEr	−3.78493 × 10−13 + 4.18932 × 10−9i	$\hphantom{-}1.71921\times 10^{-10}+7.27673\times 10^{-7}i$
δFr, x	−6.53158 × 10−10 + 2.43530 × 10−13i	−6.53158 × 10−8 + 5.44986 × 10−11i
δFr, y	$\hphantom{-}1.17047\times 10^{-7}-4.67894\times 10^{-10} i$	$\hphantom{-}1.17049\times 10^{-7}+3.10689\times 10^{-11}i$
ω	4.10391 + 2.58718 × 10−4i	4.10391 + 7.45083 × 10−4i

Download table as:  ASCIITypeset image

Table 3. Examples of Solutions for Radiation-modified Fast Mode

${\mathbb {P}}$	0.01	0.01
σ	0.01	10
δρ	10−3	10−3
δvx	$\hphantom{-}2.08441\times 10^{-3}+2.26355\times 10^{-5} i$	$\hphantom{-}1.97671\times 10^{-3}+4.44472\times 10^{-6}i$
δvy	−1.29633 × 10−3 + 3.15959 × 10−5i	−1.47025 × 10−3 + 8.22391 × 10−6i
δP	$\hphantom{-}1.64144\times 10^{-3}+1.30879\times 10^{-4} i$	$\hphantom{-}1.00104\times 10^{-3}+2.69749\times 10^{-5} i$
δBy	$\hphantom{-}2.09358\times 10^{-3}-2.82848\times 10^{-5} i$	$\hphantom{-}2.25121\times 10^{-3}-7.53016\times 10^{-6} i$
δEr	$\hphantom{-}1.89341\times 10^{-8}+7.85298\times 10^{-9} i$	$\hphantom{-}3.68258\times 10^{-6}+9.54978\times 10^{-5}i$
δFr, x	−8.33183 × 10−7 + 4.08350 × 10−6i	−1.97374 × 10−5 + 7.74322 × 10−7i
δFr, y	−1.69620 × 10−7 + 2.64650 × 10−8i	−1.96034 × 10−7 + 1.12087 × 10−9i
ω	13.0967 + 0.142223i	12.4200 + 2.79270 × 10−2i
${\mathbb {P}}$	100	100
σ	0.01	10
δρ	10−3	10−3
δvx	$\hphantom{-}1.97654\times 10^{-3}+5.75439\times 10^{-8} i$	$\hphantom{-}1.97651\times 10^{-3}+6.93247\times 10^{-6}i$
δvy	−1.47061 × 10−3 + 1.06106 × 10−7i	−1.47055 × 10−3 + 1.28333 × 10−5i
δP	$\hphantom{-}1.00000\times 10^{-3}+3.10791\times 10^{-7} i$	$\hphantom{-}9.99999\times 10^{-4}+3.17242\times 10^{-7}i$
δBy	$\hphantom{-}2.25153\times 10^{-3}-9.72688\times 10^{-8} i$	$\hphantom{-}2.25147\times 10^{-3}-1.17512\times 10^{-5} i$
δEr	−1.20333 × 10−12 + 1.26774 × 10−9i	−3.24677 × 10−9 + 1.26773 × 10−6i
δFr, x	−1.97654 × 10−9 + 1.11447 × 10−12i	−1.97803 × 10−9 + 2.44573 × 10−10i
δFr, y	−1.93101 × 10−7 + 2.39955 × 10−8i	−1.96073 × 10−7 + 1.73546 × 10−9i
ω	12.4190 + 3.61559 × 10−4i	12.4188 + 4.3558 × 10−2i
${\mathbb {P}}$	100	1
σ	0.1	10
δρ	10−3	10−3
δvx	$\hphantom{-}1.97654\times 10^{-3}+7.43855\times 10^{-8} i$	$\hphantom{-}1.97654\times 10^{-3}+2.14211\times 10^{-7}i$
δvy	−1.47061 × 10−3 + 1.37659 × 10−7i	−1.47061 × 10−3 + 3.96547 × 10−7i
δP	$\hphantom{-}1.00000\times 10^{-3}+3.42168\times 10^{-8} i$	$\hphantom{-}1.00000\times 10^{-3}+5.81556\times 10^{-7}i$
δBy	$\hphantom{-}2.25153\times 10^{-3}-1.26062\times 10^{-7} i$	$\hphantom{-}2.25153\times 10^{-3}-3.63108\times 10^{-7} i$
δEr	−1.52464 × 10−12 + 1.26774 × 10−8i	$\hphantom{-}2.54552\times 10^{-9}+2.20203\times 10^{-6} i$
δFr, x	−1.97654 × 10−9 + 2.53281 × 10−12i	−1.97654 × 10−7 + 5.86240 × 10−10i
δFr, y	−1.96050 × 10−7 + 2.45310 × 10−9i	−1.96081 × 10−7 + 7.72241 × 10−11i
ω	$\hphantom{-}12.4190+4.67378\times 10^{-4}i$	$\hphantom{-}12.419+1.34593\times 10^{-3}i$

Download table as:  ASCIITypeset image

The results of our tests for slow magnetosonic waves are shown in Figure 8, while the results for fast magnetosonic wave are shown in Figure 9. The left panel in each figure compares the numerically measured values of the phase velocity and damping rate with linear theory. Note that the solution to the dispersion relation as a function of optical depth and ${\mathbb {P}}$ for both magnetosonic modes is very similar to that for radiation-modified acoustic waves in hydrodynamics (see Figure 9). Once again, we find good agreement between our numerical values for the phase velocity versus linear theory over this parameter regime. The damping rate also shows good agreement, except when ${\mathbb {P}}=100, \tau _a <1$ and ${\mathbb {P}}=1, \tau _a>0.1$. In these cases, the damping rate is so small that our numerical measured values are not converged at this resolution. Nevertheless, as Figure 9 demonstrates, our algorithm converges to the correct solution for ${\mathbb {P}}=1, \tau _a=1$, albeit at first order. We have confirmed our algorithm also converges for the other parameters shown in this figure as well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right panels of Figures 8 and 9 show the convergence rate of the L1 error with numerical resolution for the slow and fast magnetosonic waves, respectively, in a fully 3D domain of size 2L × L × L. As before, we measure first-order convergence in a parameter regime where the solution to the radiation moment equations dominates the error, and second-order convergence in the opposite limit.

Shock tubes have long been used as a test of hydrodynamic and MHD codes in the nonlinear regime, since the structure of a planar shock can be computed analytically in these cases. However, in radiation hydrodynamics and MHD, shock structure is much more complicated, and is difficult to compute analytically. At low Mach numbers $\newmathcal {M}$, radiation can diffuse upstream of the shock, heating the gas and forming a smooth precursor in the temperature. As $\newmathcal {M}$ increases, the gas temperature near the shock front can begin to exceed the downstream value, forming a Zel'dovich spike (e.g., Zel'Dovich & Raizer 1967; Mihalas & Mihalas 1984). This spike is followed by a relaxation region where the gas temperature cools to its far downstream value. When a Zel'dovich spike is formed, if the downstream gas temperature (after the relaxation region) is larger than the value immediately upstream of the shock (at the end of the precursor region), the shock is called subcritical. On the other hand, if the downstream gas temperature after the relaxation region is the same as the upstream gas temperature, the shock is termed supercritical. Subcritical shocks are formed at lower $\newmathcal {M}$ than supercritical shocks. Finally, as $\newmathcal {M}$ is increased further, the downstream radiation pressure can exceed the gas pressure, the Zel'dovich spike disappears, and the solution becomes everywhere smooth again.

Although many authors have presented numerical solutions to radiating shocks as a test of their algorithms, the lack of analytic solutions inhibits quantitative testing. Recently Lowrie & Edwards (2008) have studied the structure of radiation-modified shocks in the non-equilibrium diffusion limit as a function of Mach number. They give a very clear explanation of how these structures can be computed by combining smooth solutions to ordinary differential equations with discontinuous jumps determined by the Rankine–Hugoniot relations when needed. These solutions not only provide interesting insight into the structure of radiating shocks, but also provide a quantitative test of time-dependent radiation hydrodynamic codes.

To test our algorithms, we follow the semi-analytic method described in (Lowrie & Edwards 2008, hereafter LE) to calculate the shock structure at a variety of Mach numbers, using a non-dimensional pre-shock solution of ρ₀ = T₀ = E_{r, 0} = 1 and $v_{0}=\newmathcal {M}$ in the upstream state. We then initialize this solution on a 1D grid by volume averaging each conserved variable to our numerical mesh. The LE solutions are given in the comoving frame, thus we transform the comoving radiation flux and energy density to the Eulerian frame, using Equation (7). In order to make our solutions match those presented in LE, we use a dimensionless speed of light ${\mathbb {C}}=1.732\times 10^3$ and pressure ${\mathbb {P}}=10^{-4}$ and absorption and scattering opacities of σ_a = σ_t = 577.4. The gas temperature T is calculated as T = P/(0.6ρ) and γ = 5/3. The Eddington approximation is used so P_r = E_r/3. We use a computational domain of size L, where L is large enough to capture the upstream precursor and downstream radiative relaxation regions (the size of these regions are given by the LE solution itself), with a grid of 1024 cells in all cases (except $\newmathcal {M}=50$ where we use 2048 cells). Input boundary conditions at the upstream values are used on the left, and outflow boundary conditions are used on the right, with the shock propagating in the negative x-direction. We then let the code evolve this solution for several flow crossing times, $t_f = L/\newmathcal {M}$. Ideally, the solution should remain stationary on the mesh. A small shift in the position of the shock front can be expected due to truncation error in the averaging of the initial solution to the hydrodynamic grid. A serious error in the algorithm or its implementation would be revealed if the code cannot hold the input solution.

We begin by presenting our numerical solution for $\newmathcal {M}=1.05$ in Figure 10. In this case there is no discontinuous jump in any variable, but rather only a smooth precursor. Next, Figure 11 shows the solution for $\newmathcal {M}=3$. Now a discontinuous shock jump has appeared, along with a Zel'dovich spike. For these parameters, this shock is subcritical. The inset panel in the plot of gas temperature shows a blow-up of the region near the spike. The numerical solution shows agreement with the semi-analytic solution of LE. Next, Figure 12 shows the solution for $\newmathcal {M}=5$. Now there is a discontinuous jump in some variables (such as density and velocity), but the gas temperature is continuous apart from the Zel'dovich spike. This solution corresponds to a supercritical shock. Again, the blow-up of the gas temperature near the spike region reveals agreement between the numerical and semi-analytic solutions. The position of the spike has moved a few cell widths, but has quickly relaxed to a steady structure at the correct shock speed. Finally, Figure 13 shows the solution for $\newmathcal {M}=50$. In this case, the radiation pressure is about 10 times gas pressure in the downstream flow. Now all variables are smooth, with no jumps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

LE have shown that the gas temperature profile near the shock front is very sensitive to small changes in $\newmathcal {M}$ when $\newmathcal {M}$ is small. To explore whether our algorithms can capture these subtle changes, Figure 14 plots the gas temperature for six Mach numbers between $\newmathcal {M}=1.05$ and 5. The values are chosen to correspond to Figures 4–11 of LE, where these changes are noted and discussed. This region covers the formation of a Zel'dovich spike, and the increase in the amplitude of the spike with increasing $\newmathcal {M}$. It is clear from the figure, in which the numerical solution is compared to the corresponding solution from LE, that our algorithm captures each phase accurately.

Zoom In Zoom Out Reset image size

Figure 14 shows that our algorithm captures the transition from subcritical to supercritical shocks accurately at low $\newmathcal {M}$, while Figure 13 shows the solution with $\newmathcal {M}=50$ (in the strongly radiation pressure dominated regime) is also captured accurately. However, we have found that solutions at $\newmathcal {M} \approx 30$ are inaccurate at the resolutions we use. As the Mach number increases beyond $\newmathcal {M}=5$, the width if the Zel'dovich spike compared to the width of the precursor region gets smaller, until the spike disappears in the radiation-pressure-dominated case near $\newmathcal {M}=50$. We find that unless the grid resolution is smaller than the thickness of the Zel'dovich spike, our algorithm will not hold a steady solution. For example, at $\newmathcal {M}=30$ this would require a resolution of about 10⁶ cells for a uniform grid. Clearly, in this case either static or adaptive mesh refinement would be useful to resolve the spike. With enough resolution to resolve the spike, our algorithms provide an accurate solution in this regime.

Calculating the structure of radiation-modified shocks without invoking the Eddington approximation or other simplifying assumptions is extremely challenging. Sincell et al. (1999) have used a time-dependent radiation hydrodynamics code to compute the structure of radiation-modified shocks. We have found that comparison of our solutions in this regime to the results reported by Sincell et al. (1999) is a useful code test.

In many ways, the numerical algorithm used by Sincell et al. (1999) is similar to that described here. They solve the radiation moment equations using a variable Eddington factor computed from a formal solution of the transfer equation, assuming LTE and gray opacities. However, the underlying hydrodynamic solver they use is quite different from that adopted in this work. Sincell et al. (1999) use methods based on artificial viscosity for shock capturing, and solve the internal rather than total gas energy equation. The Godunov methods used here use Riemann solvers rather than artificial viscosity for shock capturing, and are based on the total gas energy equation.

Since the Mach number $\newmathcal {M}$ is not known independently of the shock structure, it is not possible to compute the solution in the frame of the shock. Instead, we initialize a flow that generates a shock, and allow the shock to propagate a large enough distance to settle into a steady structure. To generate the shock, we collide two symmetric flows at the center of the domain, producing two identical shocks that propagate symmetrically away from the center. This avoids having to implement reflecting boundary conditions in the RT solver that is used to compute the VET (Davis et al. 2012) necessary if the shock is generated by reflection of the flow off a wall. The initial conditions for the flow are chosen to match those used by Sincell et al. (1999, see also SS10) for the subcritical shock solution discussed in Section 2.1 of their paper. We use a computational domain that spans −2 ⩽ x ⩽ 2. Initially, the density is one, the gas pressure is P = 0.923ρT, and the gas temperature is T = 1 + 7.5|x|/2. For x < 0, the flow velocity is v = 20, while for x > 0 it is v = −20. The dimensionless pressure and speed of light are ${\mathbb {P}}=1.08\times 10^{-10}$ and ${\mathbb {C}}=10^6$, and γ = 5/3. The radiation temperature is the same as gas temperature, the radiation flux is ${\mbox{\boldmath $F$}}_r=-(\partial E_r/\partial x)/(3\sigma _a)$ (the Eddington approximation is assumed in the initial condition only because the density is uniform initially), where the absorption opacity σ_a = 10. The VET is initialized with only non-zero diagonal components of 1/3; thereafter the VET is computed self-consistently with the flow using short characteristics.

Profiles of various quantities at t = 0.03 are shown in Figure 15 in the region 0 ⩽ x ⩽ 0.5. (We have checked, as another test of our code, that the solution is exactly symmetric with respect to x = 0.) Profiles of each variable show the characteristic structure of a subcritical shock, including a strong radiative precursor and Zel'dovich spike. Perhaps of greatest interest is the profile of the xx component of the VET shown in the lower right panel. Upstream of the shock f_xx is much bigger than 1/3, as expected. Downstream of the shock, it is close to 1/3. However, in the vicinity of the shock itself, f_xx is smaller than 1/3, in agreement with the results of Sincell et al. (1999, see their Figure 8). The reason for the behavior is discussed in Section 2.3 of Sincell et al. (1999). Most of the upstream radiation is generated within one optical depth of the shock front. Rays parallel to the shock front go through a longer path of source than rays perpendicular to the shock front. Thus, the intensity along rays parallel to the shock front is larger than those perpendicular, leading to a value of f_xx less than 1/3. Note that since the empirical relations for most flux limiters (e.g., Levermore & Pomraning 1981) bound the Eddington factor between 1/3 and one along the direction of radiation energy density gradient, FLD cannot get the right answer in this case.

Zoom In Zoom Out Reset image size

We find that the most quantitative comparison to the previously published results of Sincell et al. (1999) is to plot the gas temperature, pressure, and radiative flux against the inverse compression ratio η ≡ ρ₀/ρ, where ρ₀ is upstream density. Figure 16 shows our result at t = 0.3, for direct comparison to Figures 2 and 3 in Sincell et al. (1999). An approximate analytic expression for these quantities by Zel'Dovich & Raizer (1967) (and also given in Section 1.3 of Sincell et al. 1999) is shown as solid lines. The numerical solution of Sincell et al. (1999) is in fact in poor agreement with the analytic predictions, and these authors argue this is due to approximations made in the latter. It is clear from Figure 16 that our solutions agree very well with the approximate analytic solution, and therefore do not agree with the results of Sincell et al. (1999) for these quantities. Recall these authors use a numerical method based on artificial viscosity for shock capturing, and the internal rather than total energy equation. The plot measures the internal structure of the shock, which could be quite sensitive to both the form of viscosity used, and the degree to which energy is conserved. Our solution may be in better agreement with the analytic expectations since our algorithm uses a Riemann solver rather than artificial viscosity to capture shocks, and is based on the total energy equation.

Zoom In Zoom Out Reset image size

Photon bubble instability is an overstable mode present in radiation supported, magnetized atmospheres. It was first noticed in numerical simulations of neutron star polar cap accretion flows (Klein & Arons 1989), and was subsequently confirmed in a linear stability analysis for neutron star atmospheres (Arons 1992). More recently, Gammie (1998) and Blaes & Socrates (2003) recognized a similar instability occurs in accretion disk atmospheres, and analyzed the properties in the linear regime. Numerical simulations that investigate the nonlinear regime were performed by Turner et al. (2005) using the FLD module in the ZEUS code (Turner & Stone 2001). Since both the linear and nonlinear regimes of this instability have been studied in detail in the literature, and since fundamentally it relies on radiation-modified MHD waves, it provides an excellent test of our radiation MHD algorithms in multidimensions. In this section, we repeat one set of parameters reported by Turner et al. (2005) and compare our results with solutions from ZEUS. In order to make direct the comparisons, we adopt the Eddington approximation, so that the diagonal components of the VET are fixed to be 1/3.

The initial conditions consist of a stratified atmosphere in mechanical and thermal equilibrium, in which a uniform radiation flux balances gravity. The structure of the atmosphere is computed from the solution of the hydrostatic equilibrium equations

$$\begin{eqnarray} \frac{dP}{dy}&=&{-}\rho g+\sigma F_r,\ \ P=\rho T,\nonumber \\ \frac{1}{3}\frac{dE_r}{dy}&=&{-}\sigma F_r,\ \ E_r=T^4, \end{eqnarray} \tag{ 27 }$$

with a constant gravitational acceleration g = 1 in the $-\hat{y}$ direction, and an opacity of σ = 137.11ρ + 0.15ρ²T^−3.5 (which includes both electron scattering opacity and Kramer's opacity). The comoving vertical radiation flux $\hat{F}_{r,0}=2.52\times 10^{-4}$. These equations are integrated in a domain of size (0, 2.75) × (26.11, 28.86) in our dimensionless units, starting from the midplane y = 27.485 where the dimensionless pressure, temperature, density, and radiation energy density are all set to 1. In this test, the dimensionless parameter ${\mathbb {P}}=762$ (so the ratio between radiation pressure and gas pressure in the midplane is 254), and the dimensionless speed of light is ${\mathbb {C}}=7557.74$. The initial conditions include a uniform, horizontal magnetic field in the $+\hat{x}$ direction with magnetic pressure 10% of midplane radiation pressure. Our simulations use a resolution of 256 × 256 grid points. The instability is seeded with random perturbations of the density with amplitude 10⁻³ at the midplane. The initial condition for this simulation is similar (but not identical; the amplitude of perturbation is 10⁻⁸ in Turner et al. 2005) to the initial condition used for Figure 9 of Turner et al. (2005). This model is appropriate for the surface layers of an accretion disk at 20 Schwarzschild radii from a 10⁸ M_☉ black hole. In this case, the middle of the simulation domain corresponds to a location one scale height above the disk midplane, and the simulation domain extends from 0.95H to 1.05H. The total optical across the domain in the vertical direction is 364, thus the Eddington approximation should be an excellent description. In our dimensionless units, one orbital period is 6.28 time units.

To measure the linear growth rate of the photon bubble instability as a function of the vertical and horizontal wavenumbers k_x and k_y, we follow the method used by Turner et al. (2005). We Fourier transform the horizontal velocity v_x at 50 snapshots between times of t = 0.8 and 2.8 equally spaced at a time interval of 0.04. We assume exponential growth at each wavenumber independently, and therefore calculate the slope of the logarithm of the power with time between each neighboring time, and average the resulting rates from all the snapshots together. This approach only approximates the true growth rates in a stratified atmosphere. The result is shown in Figure 17, which can be compared directly to Figure 9 of Turner et al. (2005) (which also shows solution of linear analysis from Blaes & Socrates 2003). There is significant noise in the measurement of the growth rates, especially at small k. However, clearly the correct trends predicted by the linear analysis are reproduced, with the fastest growth occurring for modes with k_x = k_y. The noise in Figure 9 of Turner et al. (2005) is much smaller because they can use an initial amplitude as small as 10⁻⁸. This cannot be achieved in Athena because the pressure gradient and gravity cannot be balanced to round-off error (without special modification to the reconstruction algorithm) and there is always noise in the background state.

Zoom In Zoom Out Reset image size

In the nonlinear regime, Turner et al. (2005) found the instability produces shock trains that propagate from the bottom to the top of the atmosphere. With a purely horizontal background field, there is no preference for shocks propagating in either direction, resulting in a symmetric pattern of shocks in the atmosphere. The left panel of Figure 18 shows the density and velocity distribution at time 9.83 when the shock train pattern becomes relatively strong. The velocity vectors are nearly horizontal because flow is confined to be parallel to the strong horizontal magnetic field. This figure can be compared with the first panel of Figure 19 of Turner et al. (2005): our results are very similar.

Zoom In Zoom Out Reset image size

In the FLD approximation used by Turner et al. (2005), the comoving radiation flux is always assumed to be along the direction of the gradient of the comoving radiation energy density. Although we adopt the Eddington approximation in our simulation, we make no assumption about the direction of the flux. Thus, it is possible for us to test whether the radiation flux is indeed everywhere parallel to the gradient of E_r in this highly nonlinear flow. Figure 18 shows the distribution of cos  θ at t = 9.83, where the angle θ is the angle between the perturbed comoving radiation flux in our simulation ${\mbox{\boldmath $F$}}_{r,0} - \hat{F}_{r,0}$ and the direction of the flux predicted by FLD $-{\mbox{\boldmath $\nabla $}}E_{r,0}/(3\sigma _t)$. If the FLD approximation is valid, cos  θ will be 1 everywhere. We find that near the shock fronts, where there is a large density gradient, this assumption clearly fails. Since the optical depth across the domain is so large, the Eddington approximation should apply, and the flux limiter used by Turner et al. (2005) should have a numerical value very close to the value of 1/3 adopted here. Thus, the discrepancy in the direction of the radiation flux is not due to the value of the Eddington factor, but rather because in FLD the radiation flux is always along the gradient of radiation energy density.

Since the direction of the flux is not known independently of the energy density, it is well known that methods based on the diffusion approximation cannot represent shadows. This is often demonstrated through the irradiation of very optically thick structures with a beamed radiation field (e.g., Hayes & Norman 2003; González et al. 2007). In this section, we demonstrate our methods capture shadows accurately, and moreover we present the dynamical evolution of a cloud ablated by an intense radiation field. VET is crucial to get the solution correctly, which is calculated from the transfer module using short characteristic. The tests described in Davis et al. (2012, see Figure 6 of that paper) show the accuracy of the radiation transfer solver for two radiation beams. Here we show that when coupled to the VET method, our solver captures shadows correctly.

Our test is performed in a rectangular box of size (− 0.5, 0.5) × (− 0.3, 0.3) cm. Initially the background medium has density ρ₀ = 1 g cm⁻³ and temperature T₀ = 290 K. An overdense clump is located in an elliptical region r ≡ x²/a² + y²/b² ⩽ 1, with a = 0.1 cm and b = 0.06 cm. The density inside this region is ρ(x, y) = ρ₀ + (ρ₁ − ρ₀)/[1.0 + exp (10(r − 1))] with ρ₁ = 10ρ₀. The clump is in pressure equilibrium with its surroundings, so the interior is much colder than the ambient medium. The initial radiation temperature is the same as the gas temperature everywhere. Pure absorption opacity is used with σ = (T/T₀)^−3.5(ρ/ρ₀)² cm⁻¹. The radiation flux ${\mbox{\boldmath $F$}}_r$ is zero everywhere initially. Outflow boundary conditions are used on both y boundaries, and on the right x boundary. A constant radiation field with temperature T_r = 6T₀ is input through the left x boundary at angles of ±14° with respect to the x-axis. This radiation field is also input along the top y boundary at −14°, and along the bottom y boundary at +14°, in order to mimic the radiation field from an infinite plane along the left x boundary at x = −0.5. At the left x boundary, the gas temperature and density are fixed to T₀ and ρ₀, respectively. A resolution of 512 × 256 cells is used for this test. The dimensionless speed of light ${\mathbb {C}}=1.9\times 10^5$, and the parameter ${\mathbb {P}}=2.2\times 10^{-15}$.

In the ambient medium, the photon mean free path is the length of simulation domain, so that the radiation can propagate freely across the box. However, inside the clump the photon mean free path is only 3.2 × 10⁻⁶ of the length of simulation box due to the high density and low temperature. Thus, radiation coming from the left boundary cannot penetrate the clump, and shadows are cast toward the right boundary. By using an input radiation field at two angles, both umbra and penumbra are formed. This makes the test more difficult because ad hoc closures that capture only one direction for the flux will not represent both the umbra and penumbra correctly.

In Figure 19, we show the radiation energy density, and xx and yy components of the VET after one time step at δt = 1.0e⁻³. This time step corresponds to the CFL stability condition set by the adiabatic sound speed in the ambient gas and the grid resolution. The ratio of this time step to the light crossing time of the domain is about 500, thus after one time step both radiation beams have crossed the domain and a steady-state shadow structure is formed.

Zoom In Zoom Out Reset image size

Since the radiation temperature is low in this first test, the clump undergoes negligible dynamical evolution. To demonstrate the dynamics of the ablation of a clump, we have repeated this test with an increase the background temperature T₀ by a factor of 357 such that the dimensionless ${\mathbb {P}}=10^{-7}$. Then the radiation pressure of the applied radiation field at the left x boundary is about 10⁻³ of the gas pressure. This non-negligible radiation field pushes and compresses the clump from the left and the shape of the clump is changed. A shock is driven into the clump during this process. The clump is also heated up by the radiation field and outflow from the surface of the clump is formed on the left side. The density distribution and velocity field at three different times are shown in Figure 20. The adiabatic sound speed in the ambient medium is 1.3 in our units and 0.39 inside the clump. In the first panel of Figure 20, the density on the left side of the clump is increased by a factor of two due to the radiation pressure. Outflow is generated from the surface of the clump as shown by the velocity vectors. The flow is subsonic in the ambient medium and supersonic with respect to the sound speed inside the clump. In the second panel, we clearly see a shock driven to the clump with Mach number up to 3. In the third panel, the clump is completely destroyed by the radiation field and several shock fronts are formed. The ablation of clouds by strong radiation field is in fact of great interest in the interstellar medium (ISM; e.g., Bertoldi & McKee 1990; Bally 1995) and near the central region of quasars (e.g., Mathews 1983). Further explore of this flow is warranted, but beyond the scope of this paper.

Zoom In Zoom Out Reset image size

A central ingredient of the algorithm presented in this paper is that it uses a VET to close the radiation moment equations, rather than relying on the diffusion approximation. It is therefore worthwhile to discuss specific examples where there will be differences between these approaches.

The first difference is that by solving the moment equations for the radiation flux, we are able to solve for the direction of the flux self-consistently with the flow, rather than assuming it lies in the direction of ∇E_r. Both the photon bubble instability (see Figure 18) and the shadow test (see Figure 19) demonstrate the importance of keeping the direction of the flux.

A second difference is that the value of the flux limiter is always an approximation (of unknown reliability) to the true value of the diagonal components of the VET. In fact, the planar subcritical shock test (see Figure 15) shows that near the shock front, the VET can be less than 1/3 along the direction of the radiation energy density gradient, which clearly cannot be represented by the formulae used for most limiters (e.g., Levermore & Pomraning 1981). Based on the density distribution from a snapshot of simulations for black hole accretion disk done by Hirose et al. (2006), Davis et al. (2012) compare the Eddington tensor used in the FLD approximation and the VET calculated based on our transfer module. They also show significant difference (see Figure 7 of that paper). In the photosphere, the Eddington tensor given by the diffusion approximation can be much larger than the values given by the VET.

A third difference is that we do not drop the time-derivative of the flux ∂F_r/∂t in the moment equations for the flux. This term is important when the inertia of the radiation field is important, and dropping this term can lead to incorrect behavior in this regime. For example, consider a uniform medium with density ρ = 1, absorption opacity σ_a = 10, and scattering opacity σ_s = 10 in thermal equilibrium (E_r = T⁴ = 1 in our dimensionless units) with zero fixed frame radiation flux, and an initially uniform velocity along the x-direction v₀ = 1. Adopt values for the dimensionless speed of light ${\mathbb {C}}=100$ and parameter ${\mathbb {P}}=1000$, and use the Eddington approximation. In the diffusion approximation, because the radiation energy density is always uniform, the comoving radiation flux will be always be zero and the system will always be in thermal equilibrium. However, this is not true if the $\partial {\mbox{\boldmath $F$}}_r/\partial t$ term is kept. Since the gas has a non-zero velocity with respect to the radiation field, the comoving radiation flux is not zero, and there will be a drag force from the radiation field on the fluid. (a consequence of the fact that solutions to the equations of radiation hydrodynamic are not Galilean invariant). Because of the drag with the radiation field, momentum will be exchanged and the gas will be heated. In this example, this process can be described by

$$\begin{eqnarray} \frac{\partial (\rho {\mbox{\boldmath $v$}})}{\partial t} &=&\mathbb {P}\sigma _t\left({\mbox{\boldmath $F$}}_r-\frac{{\mbox{\boldmath $v$}}E_r+{\mbox{\boldmath $v$}}\cdot {\mbox{\boldmath $P$}}_r}{\mathbb {C}}\right),\ \nonumber \\ \frac{\partial {\mbox{\boldmath $F$}}_r}{\partial t}&=&-\mathbb {C}\sigma _t\left({\mbox{\boldmath $F$}}_r-\frac{{\mbox{\boldmath $v$}}E_r+{\mbox{\boldmath $v$}}\cdot {\mbox{\boldmath $P$}}_r}{\mathbb {C}}\right). \end{eqnarray} \tag{ 28 }$$

In solutions to these equations, ρ and E_r are unchanged while the change of velocity with time t can be described by

$$\begin{eqnarray} v(t)&=&v_0\left(1-\frac{{\mathbb {C}}}{{\mathbb {C}}+4{\mathbb {P}}/(3{\mathbb {C}})}\right)\exp \left[-{\mathbb {P}}\sigma _t\left(\frac{{\mathbb {C}}}{{\mathbb {P}}}+\frac{4}{3{\mathbb {C}}}\right)t\right]\nonumber\\ &&+\frac{v_0{\mathbb {C}}}{{\mathbb {C}}+4{\mathbb {P}}/(3{\mathbb {C}})}. \end{eqnarray} \tag{ 29 }$$

We compare this analytic solution for this problem, and a numerical solution generated by our algorithms, in Figure 21. Our solutions capture this momentum exchange process accurately, which is not possible with FLD. The importance of this process is characterized by the parameter ${\mathbb {P}}/{\mathbb {C}}^2$, which measures the importance of momentum carried by the radiation field compared with the material momentum.

Zoom In Zoom Out Reset image size

A fourth difference between our VET algorithm and FLD is that we keep the non-diagonal components of the VET. In the FLD approximation, the radiation pressure is always a diagonal tensor. Thus, radiation viscosity (e.g., Mihalas & Mihalas 1984; Castor 2004) cannot be captured in this approximation. However, with the VET approach, the off-diagonal components are computed self-consistently with the flow. Thus, the effects of radiation viscosity are captured naturally in simulations of accretion flows.

The importance of the differences between VET and FLD depends on the particular application of interest. For example, in the inner regions of accretion disks around supermassive black holes, we expect the momentum carried by the radiation field to be a significant fraction of material momentum (${\mathbb {P}}/{\mathbb {C}}^{2}$ is not negligible), and photon viscosity may be significant (e.g., Agol & Krolik 1998). Thus, adopting algorithms based on VET is important for our applications.

For solving the application problems of interest, our algorithms must not only be accurate, but also efficient. In this subsection, we report preliminary performance and scaling of the method on parallel systems. We continue to work on optimization of the algorithms, and it is likely substantial improvements in performance are still possible.

For reference, it is useful to compare to the performance of Athena in ideal MHD without radiation. On a single core of a 2.5 GHz Intel Xeon processor, Athena updates roughly 1 × 10⁵ grid cells per second for 3D MHD calculations with a 32³ grid. Weak scaling is very good, with an efficiency (defined as the performance per core in a parallel calculation, compared to the performance of a serial calculation) of 0.9 on up to 10⁵ processors.

With radiation hydrodynamics and MHD, there are two potential bottlenecks. The first is the implicit solution of the linear radiation subsystem, which requires inverting one matrix per time step (this is already an advantage over methods that require implicit solution of nonlinear systems which generally require many matrix inversions per time step). To invert the matrix, we have explored the use of algorithms implemented in both the LIS (e.g., Nishida 2010) and Hypre (e.g., Falgout et al. 2006) libraries, as well as our own implementation of the multigrid (e.g., Press et al. 1992) method. Both libraries provide a variety of iterative algorithms and preconditioners for solving sparse matrices, and can be used on both serial and parallel systems. The convergence rate of the matrix solver is limited by several things, such as the opacity and error tolerance.

To check the serial performance of each of these methods, we use the linear hydrodynamic wave problem as described in Section 6.1 on a 3D domain and a 32³ grid. With the LIS library on a single 2.5 GHz Intel Xeon processor with tolerance level 10⁻⁸, we achieve roughly 1.7 × 10⁴ cell updates per cpu second in optical thick regimes with the ILU preconditioner. The Hypre library is somewhat faster, achieving about 3 × 10⁴ cell updates per cpu with the GMRES solver and the BoomerAMG preconditioner. Because the sparsity pattern of the matrix we invert is not common (see Appendix A), we must use the general matrix solvers in these libraries, which have very poor cache performance. By writing our own multigrid algorithm for inverting the matrix, we can take advantage of the specific structure of our matrix to substantially improve performance. With multigrid, we can achieve 4.3 × 10⁴ cell updates per cpu for this test, and the convergence rate is independent of whether the linear waves are optically thick or thin. With both the LIS and Hypre libraries, performance is about 2–5 times slower in the optically thin regime. Currently, the multigrid solver is our workhorse algorithm.

A second bottleneck to performance in simulations with radiation is the formal solution of the transfer equation. The cost of this step depends on the number of rays required to resolve the angular dependence of the radiation field, and whether the LTE approximation can be used to compute the source function (see Davis et al. 2012). In our tests, we have found sufficient accuracy is obtained with roughly 24 angles in 3D, although this number may be highly application dependent. If absorption opacity is dominant and the LTE approximation can be used, this step consumes a negligible fraction of the overall execution time. Moreover, since the cost to compute the VET scales linearly with the number of angles, using many more angles (for example, 80 instead of 24) does not affect the overall cost per time step significantly. For example, using 80 angles increases the cost per time step by less than factor of two compared to 24 angles.

However, if scattering opacity is dominant and the non-LTE integrator must be used to compute the VET, then this step can dominate the cost. With 24 angles in a non-LTE calculation, the overall performance of the code can be decreased by a factor of up to five. Moreover, using more angles decreases the convergence rate of the iteration scheme used for non-LTE problems, thus roughly tripling the number of angles (from 24 to 80) will increase the cost per time step by more than a factor of three. There are several simple steps that can be used to increase performance in this regime, however. For example, in many cases it is possible to compute the VET only every several MHD time steps, since the geometry of the radiation field (and therefore the VET) changes only slowly compared to the flow structure.

It is interesting to compare the performance of the VET algorithm to FLD. Since we have not implemented FLD directly in Athena, we cannot perform a head-to-head comparison. However, we note that because the implicit differencing most often used with FLD (e.g., Turner & Stone 2001; Zhang et al. 2011) involves nonlinear terms, Newton–Raphson iterations are required with one matrix solver per iteration. Depending on the convergence rate of the iterations, this could require several to many matrix solvers per time step. Thus, Turner & Stone (2001) reported that in the optically thin limit the FLD module in ZEUS made the code run 10 × slower, whereas we find that the VET algorithm (including the short-characteristics RT solver) is only 2–3 × slower than the standard MHD integrators in Athena. ZEUS itself is several times faster than Athena, so the VET may not in fact be faster than FLD, but it is at least comparable in performance. Using FLD because it is faster does not seem to be supported by these comparisons.

The parallel performance of the code is also limited by the matrix solver we use. Generally, in parallel calculations, the iterative linear system solver requires more steps to converge compared to the serial case. Using the linear hydrodynamic wave problem as described above, with 32³ grids on each processor, we have performed a weak scaling test for both the LIS and Hypre libraries, and our own multigrid solver. With 8 cores, the efficiency of the two libraries are about the same: 0.65, and on 256 cores this drops to 0.25. With our multigrid solver, the efficiency on 8 cores is 0.76, and drops to 0.56 on 256 cores. While this is substantially lower than the efficiency of the MHD integrator alone, it is still sufficient to enable large 3D applications. However, clearly further effort to improve performance is warranted.