MOCMC: Method of Characteristics Moment Closure, a Numerical Method for Covariant Radiation Magnetohydrodynamics

The governing equations of covariant MHDs, written in conservative form in a coordinate basis, are (Anile 1989; Komissarov 1999; Gammie et al. 2003) conservation of particle number

$$\begin{eqnarray}&&{\left(\sqrt{-g}\rho {u}^{t}\right)}_{,t}=-{\left(\sqrt{-g}\rho {u}^{i}\right)}_{,i},\end{eqnarray} \tag{ 3 }$$

conservation of four-momentum (note that we are deferring discussion of radiation-matter source terms to Section 2.2.1)

$$\begin{eqnarray}&&{\left(\sqrt{-g}{T}^{t}{}_{\nu }\right)}_{,t}=-{\left(\sqrt{-g}{T}^{i}{}_{\nu }\right)}_{,i}+\sqrt{-g}{T}^{\kappa }{}_{\lambda }{{\rm{\Gamma }}}^{\lambda }{}_{\nu \kappa },\end{eqnarray} \tag{ 4 }$$

conservation of magnetic flux

$$\begin{eqnarray}&&{\left(\sqrt{-g}{B}^{i}\right)}_{,t}=-{\left(\sqrt{-g}\left[{b}^{j}{u}^{i}-{b}^{i}{u}^{j}\right]\right)}_{,j},\end{eqnarray} \tag{ 5 }$$

and the no-monopoles constraint

$$\begin{eqnarray}&&{\left(\sqrt{-g}{B}^{i}\right)}_{,i}=0.\end{eqnarray} \tag{ 6 }$$

The MHD stress–energy tensor is

$$\begin{eqnarray}&&{T}^{\mu }{}_{\nu }=\left(\rho +{u}_{{\rm{g}}}+{P}_{{\rm{g}}}+{b}^{2}\right){u}^{\mu }{u}_{\nu }\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&+\,\left({P}_{{\rm{g}}}+\displaystyle \frac{{b}^{2}}{2}\right){\delta }^{\mu }{}_{\nu }-{b}^{\mu }{b}_{\nu }.\end{eqnarray} \tag{ 8 }$$

In the above, ρ is the fluid rest-mass density, u^μ is the fluid four-velocity, Bⁱ is the magnetic field three-vector, b^μ is the magnetic field four-vector, and $\sqrt{-g}$ is the determinant of the metric. These equations require an equation of state; throughout this work, we adopt ${P}_{{\rm{g}}}=(\gamma -1){u}_{{\rm{g}}}$, with ${u}_{{\rm{g}}}$ the gas internal energy density, although introducing more sophisticated equations of state in this framework is conceptually straightforward (e.g., Miller et al. 2019a).

The equation of radiation transport, written in invariant form (each quantity in parentheses is invariant), is (Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&\displaystyle \frac{D}{{ds}}\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right)=\left(\displaystyle \frac{{j}_{\nu }^{{\rm{a}}}}{{\nu }^{2}}\right)+\left(\displaystyle \frac{{j}_{\nu }^{{\rm{s}}}}{{\nu }^{2}}\right)-\left(\nu {\alpha }_{\nu }^{{\rm{a}}}\right)\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right)-\left(\nu {\alpha }_{\nu }^{{\rm{s}}}\right)\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right),\end{eqnarray} \tag{ 9 }$$

where D/ds is the convective derivative in phase space, s is the path length, ${j}_{\nu }^{{\rm{a}}}$ is the specific thermal emissivity, ${\alpha }_{\nu }^{{\rm{a}}}$ is the specific absorption coefficient, ${j}_{\nu }^{{\rm{s}}}$ is the effective specific emission coefficient due to scattering (Mihalas & Mihalas 1984), and ${\alpha }_{\nu }^{{\rm{s}}}$ is the specific absorption coefficient due to scattering. The specific intensity I_ν is related to the distribution function f as ${I}_{\nu }={h}^{4}{\nu }^{3}f/{c}^{2}$. We will return to source terms in Section 2.2.1; for now, it is sufficient to note that, absent interaction terms, Equation (9) preserves the invariant intensity I_ν/ν³ along characteristics.

Each characteristic is described by a position x^μ and a momentum p^μ defined such that −p^μu_μ = hν, where ν is the frequency measured by an observer with four-velocity u^μ. These quantities evolve according to the geodesic equation (note that the coordinate time t = x⁰):

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}^{i}}{{{dx}}^{0}}=\displaystyle \frac{{p}^{i}}{{p}^{0}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dp}}_{\mu }}{{{dx}}^{0}}=-\displaystyle \frac{1}{2{p}^{0}}\displaystyle \frac{\partial {g}^{\nu \lambda }}{\partial {x}^{\mu }}{p}_{\nu }{p}_{\lambda }\end{eqnarray} \tag{ 11 }$$

where in the first equation we have implicitly divided out the affine parameter, defined by dx^μ/dλ = p^μ.

For Monte Carlo particles, like superphotons in bhlight (Ryan et al. 2015), one would typically evolve the momentum p^μ directly. However, we will discretize I_ν/ν³ in frequency along each characteristic, so in practice we evolve a direction vector whose normalization has no special meaning. For a particular momentum, the direction vector is n^μ = p^μ/(hν). Effectively, we normalize n^μ such that n₀ = −1; the frequency ν of each bin i according to an observer moving with four-velocity u^μ is then

$$\begin{eqnarray}&&\nu =-\displaystyle \frac{{u}^{\mu }{n}_{\mu }{p}_{0}^{i}}{h},\end{eqnarray} \tag{ 12 }$$

where p₀ⁱ is an array of timelike components of covariant four-momenta—or for our purposes, frequencies at asymptotic spatial infinity common to all samples.³ This array defines the range of the frequency discretization. We discretize these frequencies logarithmically; note that ${\rm{\Delta }}\mathrm{log}\nu$ is unaffected by frame transformations. In terms of x^μ and n^μ, the geodesic equation is then

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}^{i}}{{{dx}}^{0}}=\displaystyle \frac{{n}^{i}}{{n}^{0}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{\mu }}{{{dx}}^{0}}=-\displaystyle \frac{1}{2{n}^{0}}\displaystyle \frac{\partial {g}^{\nu \lambda }}{\partial {x}^{\mu }}{n}_{\nu }{n}_{\lambda }.\end{eqnarray} \tag{ 14 }$$

The moments of the radiation field evolve due to advection (i.e., ignoring radiation-matter source terms) in a similar manner to the MHD stress–energy tensor,

$$\begin{eqnarray}&&{\left(\sqrt{-g}{R}^{t}{}_{\nu }\right)}_{,t}=-{\left(\sqrt{-g}{R}^{i}{}_{\nu }\right)}_{,i}+\sqrt{-g}{R}^{\kappa }{}_{\lambda }{{\rm{\Gamma }}}^{\lambda }{}_{\nu \kappa },\end{eqnarray} \tag{ 15 }$$

where ${R}^{\mu }{}_{\nu }$ is related to the specific intensity I_ν as

$$\begin{eqnarray}&&{R}^{\mu }{}_{\nu }=\displaystyle \frac{1}{{h}^{4}}\displaystyle \int \displaystyle \frac{{d}^{3}p}{\sqrt{-g}{p}^{t}}{p}^{\mu }{p}_{\nu }\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right).\end{eqnarray} \tag{ 16 }$$

However, unlike in the MHD case where we assumed a thermal, isotropic distribution of particles in the fluid frame, there is no general analytic expression for ${R}_{\,j}^{i}$. In the fluid frame, this uncertainty can be parameterized as (recall that parentheses indicate an orthonormal tetrad, and here the comoving frame of the fluid) ${R}^{(i)}{}_{(j)}=-{R}^{(0)}{}_{(0)}{\pi }^{(i)}{}_{(j)}$, where ${\pi }^{(i)}{}_{(j)}$ is the Eddington tensor, which is calculated from the specific intensity:

$$\begin{eqnarray}&&{\pi }^{(i)}{}_{(j)}=\displaystyle \frac{\displaystyle \int d\nu d{\rm{\Omega }}{I}_{\nu }{n}^{(i)}{n}_{(j)}}{\displaystyle \int d\nu d{\rm{\Omega }}{I}_{\nu }},\end{eqnarray} \tag{ 17 }$$

where dΩ is the differential solid angle in an orthonormal tetrad. Given ${\pi }^{(i)}{}_{(j)}$ and ${R}^{0}{}_{\nu }$, we can construct the entire radiation stress–energy tensor in any frame.

For the special case of very optically thick flows, the specific intensity goes to the Planck function ${B}_{\nu }$, and the Eddington tensor becomes

$$\begin{eqnarray}&&{\pi }^{(i)}{}_{(j)}=\displaystyle \frac{{{\delta }^{(i)}}_{(j)}}{3}.\end{eqnarray} \tag{ 18 }$$

When optical depths are not very large, ${\pi }^{(i)}{}_{(j)}$ is in general subject to few a priori constraints, and one must be able to evaluate Equation (17); that is, one must have knowledge of ${I}_{\nu }$.

2.2.1. Fluid–Radiation Interactions

Microphysical processes (emission, absorption, and scattering) lead to an exchange of four-momentum between the fluid and the radiation, along with a change in specific intensity along characteristics.

The fluid and radiation stress energy tensors communicate through exchange of four-momentum. We can rewrite the divergence of the radiation stress–energy tensor as a four-force density,

$$\begin{eqnarray}&&{G}_{\nu }\equiv -{R}^{\mu }{}_{\nu ;\mu }.\end{eqnarray} \tag{ 19 }$$

From Equation (2), we then have the interaction source terms on our fluid and radiation four-momenta:

$$\begin{eqnarray}&&{\left(\sqrt{-g}{T}^{t}{}_{\nu }\right)}_{,t}=\sqrt{-g}{G}_{\nu }\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\left(\sqrt{-g}{R}^{t}{}_{\nu }\right)}_{,t}=-\sqrt{-g}{G}_{\nu }.\end{eqnarray} \tag{ 21 }$$

In the fluid frame, the four-force density is (Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&{G}_{(\mu )}=\displaystyle \int d\nu d{\rm{\Omega }}\left({\alpha }_{\nu }^{{\rm{a}}}{I}_{\nu }+{\alpha }_{\nu }^{{\rm{s}}}{I}_{\nu }-{j}_{\nu }^{{\rm{a}}}-{j}_{\nu }^{{\rm{s}}}\right){n}_{(\mu )}.\end{eqnarray} \tag{ 22 }$$

We now rewrite this equation to produce something more amenable to stable integration. In particular, we want to use the I_ν to average source terms, rather than compute four-forces along individual characteristics as in a Monte Carlo approach, to reduce noise in the fluid. We will also introduce an angle-averaged approach to scattering.

The term involving ${j}_{\nu }^{{\rm{a}}}$ in the integrand does not depend on the intensity; we may integrate this emissivity directly (J^a). In addition, we can without approximation rewrite the absorption terms as frequency- and angle-averaged opacities multiplying the comoving radiation four-momentum:

$$\begin{eqnarray}&&{G}_{(0)}=-{\overline{\alpha }}^{{\rm{a}}}{R}^{(0)}{}_{(0)}-{\overline{\alpha }}^{{\rm{s}}}{R}^{(0)}{}_{(0)}-{J}^{{\rm{a}}}-\displaystyle \int d\nu d{\rm{\Omega }}{j}_{\nu }^{{\rm{s}}},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{G}_{(i)}={\overline{\alpha }}^{{\rm{a}}}{R}^{(0)}{}_{(i)}+{\overline{\alpha }}^{{\rm{s}}}{R}^{(0)}{}_{(i)}-\displaystyle \int d\nu d{\rm{\Omega }}{j}_{\nu }^{{\rm{s}}}{n}_{(i)},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{\overline{\alpha }}^{{\rm{a}}}\equiv \displaystyle \frac{1}{I}\displaystyle \int d\nu d{\rm{\Omega }}{\alpha }_{\nu }^{{\rm{a}}}{I}_{\nu },\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\overline{\alpha }}^{{\rm{s}}}\equiv \displaystyle \frac{1}{I}\displaystyle \int d\nu d{\rm{\Omega }}{\alpha }_{\nu }^{{\rm{s}}}{I}_{\nu },\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&I\equiv \displaystyle \int d\nu d{\rm{\Omega }}{I}_{\nu }.\end{eqnarray} \tag{ 27 }$$

We now consider the term involving ${j}_{\nu }^{{\rm{s}}}$ in Equation (23). For any (dΩ, dν), ${j}_{\nu }^{{\rm{s}}}$ can in principle depend on every other (dΩ, dν); all characteristics can scatter into each other. This is numerically very expensive unless treated in a probabilistic manner (e.g., Dolence et al. 2009; Ryan et al. 2015). Here, we wish to preserve our continuum approach. Therefore, we instead derive a specific scattering emissivity from the evolution of an angle-averaged intensity given by the Kompaneets equation.

We approximate scattering by considering the angle-integrated transfer equation in the comoving frame with only source terms due to scattering:

$$\begin{eqnarray}&&\displaystyle \frac{d{{ \mathcal I }}_{\nu }}{{ds}}={{ \mathcal J }}_{\nu }^{{\rm{s}}}-{{ \mathcal A }}_{\nu }^{{\rm{s}}}{{ \mathcal I }}_{\nu },\end{eqnarray} \tag{ 28 }$$

where script letters (${ \mathcal I }$, ${ \mathcal J }$, ${ \mathcal A }$) indicate a solid angle integral average. By evaluating $d{{ \mathcal I }}_{\nu }/{ds}$, we can solve for ${{ \mathcal J }}_{\nu }^{{\rm{s}}}$, and then approximate the comoving frame radiation four-force as

$$\begin{eqnarray}&&{G}_{(0)}=-{\overline{\alpha }}^{{\rm{a}}}{R}^{(0)}{}_{(0)}-{\overline{\alpha }}^{{\rm{s}}}{R}^{(0)}{}_{(0)}-{J}^{{\rm{a}}}-4\pi {{ \mathcal J }}^{{\rm{s}}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{G}_{(i)}={\overline{\alpha }}^{{\rm{a}}}{R}^{(0)}{}_{(i)}+{\overline{\alpha }}^{{\rm{s}}}{R}^{(0)}{}_{(i)}.\end{eqnarray} \tag{ 30 }$$

One obvious consequence of our angle-integrated procedure is that the integral over ${j}_{\nu }^{{\rm{s}}}$ does not contribute to G_(i), because ${j}_{\nu }^{{\rm{s}}}$ is isotropic in the fluid frame. The error associated with this procedure depends on the differential cross section; for an isotropic scattering process, this is exact. The approximation we make to scattering is separate from any treatment of emission and absorption; those terms are solved exactly. Our approximate approach to scattering is not a fundamental requirement of MOCMC. We only invoke this approximation here for computational expediency and because this treatment is likely sufficient for some applications of immediate interest.

Here, $d{{ \mathcal I }}_{\nu }/{ds}$ is zero for elastic scattering; in general, however, it depends on the scattering process under consideration. In order to proceed, we now specialize to electron scattering and introduce the Kompaneets equation (e.g., Rybicki & Lightman 1979):

$$\begin{eqnarray}&&\displaystyle \frac{\partial n}{\partial \tau }={n}_{{\rm{e}}}{\sigma }_{{\rm{T}}}c\left(\displaystyle \frac{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}{{m}_{{\rm{e}}}{c}^{2}}\right)\displaystyle \frac{1}{{x}^{2}}\displaystyle \frac{\partial }{\partial x}\left[{x}^{4}\left(\displaystyle \frac{{dn}}{{dx}}+n+{n}^{2}\right)\right],\end{eqnarray} \tag{ 31 }$$

where n_e is the electron number density, σ_T is the Thomson cross section, T_e is the electron temperature, m_e is the electron rest mass, $x=h\nu /({k}_{{\rm{B}}}{T}_{{\rm{e}}})$, $n={c}^{2}{{ \mathcal I }}_{\nu }/(2h{\nu }^{3})$ is the photon occupation number, and τ is here the proper time in the fluid frame. The Kompaneets equation is an angle-integrated (i.e., consistent with our scattering approximation) expansion of the Compton scattering kernel in the dimensionless energy transferred to photons per scattering event, $h{\rm{\Delta }}\nu /({k}_{{\rm{B}}}{T}_{{\rm{e}}})$. This value is small for nonrelativistic electrons, T_e ≲ 10⁸ K, and when small, Compton scattering becomes a diffusive flux in momentum space. At higher temperatures, one must take an integrodifferential approach to Compton scattering (e.g., Jones 1968; Aharonian & Atoyan 1981; Coppi & Blandford 1990; Dolence et al. 2009; Suleimanov et al. 2012; Ryan et al. 2015; Narayan et al. 2017; Kinch et al. 2019).

Up to a multiplicative constant ${ \mathcal C }$, $\partial n/\partial t\sim d{{ \mathcal I }}_{\nu }/{ds}$. We can therefore use the Kompaneets equation to evaluate the change in intensity, and solve Equation (28) for the angle-averaged scattering emissivity:

$$\begin{eqnarray}&&{{ \mathcal J }}_{\nu }^{{\rm{s}}}={ \mathcal C }\displaystyle \frac{\partial n}{\partial t}+{{ \mathcal A }}_{\nu }^{{\rm{s}}}{{ \mathcal I }}_{\nu }.\end{eqnarray} \tag{ 32 }$$

Equations (29) and (30) describe the exchange of four-momentum between fluid and radiation, but they do not update specific intensities. These are calculated by solving the transfer equation along each characteristic in invariant form. For each I_ν, we solve Equation (9), with ${j}_{\nu }^{{\rm{s}}}$ approximated by ${{ \mathcal J }}_{\nu }^{{\rm{s}}}$.

Our method for integrating the GRMHD equations is harm (Gammie et al. 2003), a flux-conservative shock capturing scheme. All gas and magnetic field variables are zone-centered. Second-order accuracy in time is achieved with a midpoint method. The primitive variables are

$$\begin{eqnarray}&&\left(\rho ,{u}_{{\rm{g}}},{\tilde{u}}^{1},{\tilde{u}}^{2},{\tilde{u}}^{3},{B}^{1},{B}^{2},{B}^{3}\right),\end{eqnarray} \tag{ 33 }$$

where Bⁱ is the magnetic field three-vector and ${\tilde{u}}^{i}$ is related to the spatial components of the four-velocity, but more amenable to variable inversion (McKinney & Gammie 2004). The conserved variables U are

$$\begin{eqnarray}&&\sqrt{-g}\left(\rho {u}^{0},{T}^{0}{}_{0}-\rho {u}^{0},{T}^{0}{}_{1},{T}^{0}{}_{2},{T}^{0}{}_{3},{B}^{1},{B}^{2},{B}^{3}\right)\end{eqnarray} \tag{ 34 }$$

and the fluxes Fⁱ are

$$\begin{eqnarray}&&\begin{array}{l}\sqrt{-g}(\rho {u}^{i},{T}_{0}^{i}-\rho {u}^{i},{T}_{1}^{i},{T}_{2}^{i},{T}_{3}^{i},\\ \quad {b}^{1}{u}^{i}-{b}^{i}{u}^{1},{b}^{2}{u}^{i}-{b}^{i}{u}^{2},{b}^{3}{u}^{i}-{b}^{i}{u}^{3}).\end{array}\end{eqnarray} \tag{ 35 }$$

Note that we have subtracted the rest mass from the time component of the stress–energy; this allows for greater accuracy in the internal energy of the fluid in certain cases.

We use monotonized central (second-order) or WENO5 (fifth-order; Liu et al. 1994; Tchekhovskoy et al. 2007) methods to reconstruct primitive variables at zone faces. These in turn are used to calculate left and right fluxes F_L and F_R and conserved variables U_L and U_R. The local Lax–Friedrichs approximate Riemann solver is then used to compute intercell fluxes of conserved variables,

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{g}}}^{i}=\displaystyle \frac{{F}_{{\rm{L}}}^{i}+{F}_{{\rm{R}}}^{i}-{c}_{\mathrm{top},{\rm{g}}}({U}_{{\rm{R}}}-{U}_{{\rm{L}}})}{2},\end{eqnarray} \tag{ 36 }$$

where ${c}_{\mathrm{top},{\rm{g}}}$ is the maximum fluid wavespeed in the coordinate frame.

The no-monopoles condition is enforced via flux-CT (Tóth 2000). This method is robust and preserves a numerical discretization of Equation (6) to machine precision, although approaches using upwinded electromotive forces can deliver superior performance on at least certain problems (Gardiner & Stone 2008; White et al. 2016) and can be extended to grids that are not logically Cartesian (Duffell 2016).

The conserved radiation four-momentum $\sqrt{-g}{R}^{0}{}_{\nu }$ is discretized spatially at zone centers. The samples in each zone will subsequently close the evolution equations for $\sqrt{-g}{R}^{0}{}_{\nu }$ by providing ${R}^{i}{}_{j}$. In this section, we assume that the complete stress tensor is known.

Similarly to the fluid evolution, we reconstruct from zone centers to faces using monotonized central or WENO5, and then calculate fluxes using the local Lax–Friedrichs solver. However, our radiation moment procedure differs from our treatment of advection for MHD in two important ways.

First, we reconstruct conserved variables and fluxes directly, rather than primitive variables. While reconstructing primitive variables is sometimes advantageous for hydrodynamic methods (preventing, for example, negative pressures), we have not found our alternative, more convenient approach to behave pathologically in the problems we have considered. In addition, transforming vector/tensor quantities like radiation fluxes or the Eddington tensor to locally orthonormal frames and then reconstructing them is not a unique process. While some choices of frames are obviously better than others from the standpoint of truncation error, free rotations in neighboring frames reconstructed to faces could introduce errors. Second, for the purposes of numerical diffusion, we assume that the wavespeed is always c. This is accurate for free streaming along one direction, but for nearly isotropic radiation, this will overestimate wavespeeds by a factor $\sqrt{3}$, leading to some additional numerical diffusion.

Our fluxes are then given by (again, using local Lax–Friedrichs):

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{r}}}^{i}=\sqrt{-g}\displaystyle \frac{{R}^{i}{}_{\nu ,{\rm{L}}}-{R}^{i}{}_{\nu ,{\rm{R}}}+c\left({R}^{0}{}_{\nu ,{\rm{R}}}-{R}^{0}{}_{\nu ,{\rm{L}}}\right)}{2}.\end{eqnarray} \tag{ 37 }$$

This approach can produce unacceptably large numerical diffusion when the optical depth per zone is very large: the diffusion applied by local Lax–Friedrichs overestimates that required for stability in what is essentially a parabolic problem, and Equation (37) fails to recover the asymptotic diffusion limit. Previous authors (Jin & Levermore 1996; Sa̧dowski et al. 2013; Foucart et al. 2015; Skinner et al. 2019) have treated this issue by calculating a separate flux valid in the optically thick limit, ${{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{i}$, and interpolating between these two based on the local optical depth. Here, we follow the procedure described in Foucart et al. (2015). Essentially, we split the energy and momentum fluxes into terms due to radiation advection and radiation diffusion, upwind the advection term, and average the diffusion term to construct ${{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{i}$, which we then smoothly interpolate toward from ${{ \mathcal F }}_{{\rm{r}}}^{i}$, based on the intensity-weighted frequency-averaged optical depth in the zone Δτ_zone.

The lab-frame energy and momentum fluxes in the Eddington closure, ${{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{i}={R}_{\mu }^{i}$, with ${u}_{{\rm{r}}}\equiv {R}_{(0)}^{(0)}$, is (Farris et al. 2008):

$$\begin{eqnarray}&&{R}_{\mu }^{i}=\displaystyle \frac{4}{3}{u}_{{\rm{r}}}{u}^{i}{u}_{\mu }+{F}^{i}{u}_{\mu }+{u}^{i}{F}_{\mu },\end{eqnarray} \tag{ 38 }$$

where F^μ is the radiation flux vector, defined such that F^μu_μ = 0. We treat the first term as the advection of radiation energy by the fluid, and the subsequent terms as the diffusion of radiation. We adopt the relativistic expression for the flux (but drop time derivatives and the four-acceleration term),

$$\begin{eqnarray}&&{F}_{i}=-\displaystyle \frac{1}{3\chi \rho }\left({\delta }_{i}^{j}+{u}_{i}{u}^{j}\right){u}_{{\rm{r}},j},\end{eqnarray} \tag{ 39 }$$

and evaluate F₀ from the condition ${F}_{\mu }{u}^{\mu }=0$. We then construct a final, stable, asymptotic-preserving flux ${{ \mathcal F }}_{{\rm{r}},\mathrm{asym}}^{i}$ by interpolating based on the optical depth in the zone,

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{r}},\mathrm{asym}}^{i}=\epsilon {{ \mathcal F }}_{{\rm{r}}}^{i}+(1-\epsilon ){{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{i},\end{eqnarray} \tag{ 40 }$$

where

$$\begin{eqnarray}&&\epsilon =\tanh \left(\displaystyle \frac{1}{{\rm{\Delta }}{\tau }_{\mathrm{zone}}}\right).\end{eqnarray} \tag{ 41 }$$

For the advective term in ${{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{i}$, we reconstruct the fluid primitive variables and u_r to each face. We use the reconstructed fluid primitive variables to evaluate fluid coordinate velocities. The average of the left and right coordinate velocities is used to define both the upwind direction and the advective velocity. The advection term is then evaluated using this average velocity and the reconstructed u_r on the upwind side of the face.

For the diffusive terms ${F}^{i}{u}_{\mu }+{u}^{i}{F}_{\mu }$ in a flux ${{ \mathcal F }}_{{\rm{r}},\,\mathrm{diff}}^{j}$, we first approximate ${u}_{{\rm{r}},i}$ along direction j as

$$\begin{eqnarray}&&{u}_{{\rm{r}},i}\approx \displaystyle \frac{{u}_{{\rm{r}}}^{j}-{u}_{{\rm{r}}}^{j-1}}{{\rm{\Delta }}{x}^{i}}.\end{eqnarray} \tag{ 42 }$$

We then evaluate the rest of the diffusive terms using left- and right-state fluid quantities, and take their average.

We discretize the invariant specific intensity I_ν/ν³ with a swarm of samples. Each sample has a unique position x^μ and direction vector n_μ, and carries an array of I_ν/ν³ discretized logarithmically in p₀. In practice, we initialize the positions and direction vectors by sampling uniformly in space and on the unit sphere in momentum space, respectively, but this is not a requirement.

For stationary spacetimes (such as Minkowski space and rotating black holes) p₀ is invariant (and is equivalent to the frequency at infinity for asymptotically flat spacetimes). This would not hold in dynamical spacetimes, such as compact object mergers. In a tetrad frame with coordinate four-velocity u^μ, the comoving frequency is ∝u^μn_μ. Due to our logarithmic discretization in k₀, when considering frequency bins, we evaluate Δν as $\nu {\rm{\Delta }}\mathrm{log}\nu$, where ${\rm{\Delta }}\mathrm{log}\nu$ is a constant.

We evolve x^μ and n_μ by directly integrating Equation (14), a set of ordinary differential equations, similarly to Dolence et al. (2009). We use the second-order-accurate Heun's method. We adaptively refine our integration to ensure some tolerance is met in Δn^μn_μ at each step; for the null geodesics we consider, n^μn_μ = 0, but this will not generally be conserved by our numerical geodesic integration. For spacetimes symmetric in x^μ, the source terms on n_μ are zero; these quantities are conserved.

Source terms are applied to characteristics in an operator-split fashion after the geodesic update to x^μ and n_μ just described; see Section 3.8.

We employ two frames in this work: the coordinate frame, and a set of orthonormal tetrad frames comoving with the fluid. Essentially, the transport operators are evaluated in the coordinate frame, and source terms and the Eddington tensor are evaluated in comoving frames. These tetrads are constructed with Gram–Schmidt orthogonalization (e.g., Dolence et al. 2009), producing transformation matrices between tetrad and coordinate frames:

$$\begin{eqnarray}&&{e}^{(\mu )}{}_{\mu },\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&e{}_{(\mu )}{}^{\mu },\end{eqnarray} \tag{ 44 }$$

accurate to roundoff error.

Fluid tetrads are constructed at the center of each zone. However, when boosting the samples inside a zone into that fluid frame, we construct a separate tetrad transformation at the spatial coordinate of that sample, and then simply interpret the axes of that tetrad as being "close enough" to the zone-centered fluid frame. This ensures that samples remain normalized in the zone-centered fluid frame; transforming a vector with a tetrad transformation evaluated at a different spatial coordinate will in general affect the normalization of that vector. However, this process is a source of error in sample positions on the unit sphere in the comoving frame. Note that, in curved space, there is no unique way to assign angles between vectors that do not share spatial coordinates.

In the comoving frame, we wish to compute integrals of the specific intensity over solid angle. We treat the samples as support points and then define a quadrature rule. A simple approach would be the traditional Monte Carlo method, in which every particle has equal angular "weight" and the sample intensities are simply summed over. Here, we adopt a more accurate approach, although MOCMC is largely agnostic in this respect.

To reconstruct momentum space, we construct a Delaunay triangulation, with the samples in each zone acting as vertices. To do this, we actually construct the convex hull, the smallest-volume region that contains the samples (or the closed surface one gets from "shrinkwrapping" the samples) in $\sim { \mathcal O }({N}_{\mathrm{samp}}\mathrm{log}{N}_{\mathrm{samp}})$ time using the CGAL⁴ library's implementation (Hert & Schirra 2018) of the quickhull algorithm (Barber et al. 1996). Here, N_samp is the number of samples in a zone. The facets of this hull, projected onto the unit sphere, are equivalent to the spherical Delaunay triangulation. Figure 1 shows the convex hull of 64 points sampled uniformly on the unit sphere.

Zoom In Zoom Out Reset image size

For each spherical triangle, we calculate the spherical excess e, or solid angle ${\rm{\Delta }}{\rm{\Omega }}$ subtended by the spherical triangle, with l'Huilier's theorem,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\rm{\Omega }} & = & 4{\tan }^{-1}\\ & & \times \left[\sqrt{\tan \left(\displaystyle \frac{s}{2}\right)\tan \left(\displaystyle \frac{s-\alpha }{2}\right)\tan \left(\displaystyle \frac{s-\beta }{2}\right)\tan \left(\displaystyle \frac{s-\gamma }{2}\right)}\right],\end{array}\,\end{eqnarray} \tag{ 45 }$$

where $s=\left(\alpha +\beta +\gamma \right)/2$ and α, β, γ are the angles between each pair of vertices measured from the origin. We will use these ΔΩ both to calculate the Eddington tensor and to angle-average intensities and opacities when evaluating fluid-radiation interactions.

For computational efficiency when evaluating angular quadratures over this triangulation, we interpolate the comoving-frame specific intensities for each sample onto a common frequency grid.

In the event that a zone contains fewer than four samples, the convex hull is not defined; instead, we set ${\pi }^{(i)}{}_{(j)}={\delta }^{(i)}{}_{(j)}/3$ and ${\rm{\Delta }}{\rm{\Omega }}=4\pi /{N}_{\mathrm{samp}}$. In practice, this is rare; our adaptive approach to sampling (Section 3.7) generally prevents this for even a modest (∼16) average number of samples per zone.

The evaluation of ${R}^{i}{}_{j}$ is divided into two parts. First, we calculate the comoving frame ${\pi }^{(i)}{}_{(j)}$ by integrating the frequency-integrated intensities at the vertices of the triangulation over the unit sphere. Second, we calculate the comoving frame four-momentum such that, when we transform the comoving frame radiation stress–energy tensor back to the coordinate frame, we recover the coordinate frame four-momentum we started with. Note that the evaluation of ${R}^{(0)}{}_{(\mu )}$ from ${R}^{0}{}_{\mu }$ depends on ${\pi }^{(i)}{}_{(j)}$, hence the complexity of this second part.

For integrating the Eddington tensor, we adopt a simple quadrature rule. For a spherical triangle with solid angle ΔΩ and vertices with weights (here, frequency-integrated intensities) w¹, w², and w³, the contribution to a component of the Eddington tensor is

$$\begin{eqnarray}&&{\rm{\Delta }}{\pi }^{(i)}{}_{(j)}={\rm{\Delta }}{\rm{\Omega }}\displaystyle \frac{{\sum }_{m}^{3}\tfrac{{w}_{m}}{3}{n}^{(i)}{n}_{(j)}}{{\sum }_{m}^{3}\tfrac{{w}_{m}}{3}},\end{eqnarray} \tag{ 46 }$$

where the n⁽ⁱ⁾, n_(j) are evaluated at each vertex in the sum.

To increase accuracy at modest additional cost, we adopt the method of Boal & Sayas (2004), albeit with a different quadrature rule for individual triangles. For each initial triangular face, one may subdivide this triangle into four subfaces and integrate these subfaces individually. The resulting integral of a quantity over the sphere at a refinement level N is then denoted I_N. Boal & Sayas (2004) conjecture that these refinement levels may then be combined in a Richardson extrapolation-like manner to produce higher-order results, at least for weights known exactly (we simply average our weights to midpoints, although the n⁽ⁱ⁾, n_(j) we use are exact). We use one refinement level, i.e., the coarsest realization of this approach, which is conjectured to behave as:

$$\begin{eqnarray}&&{\displaystyle \int }_{{\rm{\Omega }}}=\,\displaystyle \frac{4}{3}{I}_{N+1}-\displaystyle \frac{1}{3}{I}_{N}+{ \mathcal O }\left(\langle {\rm{\Delta }}{\rm{\Omega }}{\rangle }^{4}\right),\end{eqnarray} \tag{ 47 }$$

which is consistent with our limited numerical experiments. On nearly isotropic multizone problems in MOCMC, we have observed reductions in error in gas temperature of up to 30× compared to results based on the single level integration scheme.

Given a set of triangles, the Eddington tensor is then

$$\begin{eqnarray}&&{\pi }^{(i)}{}_{(j)}=\displaystyle \frac{{\sum }_{l}{\rm{\Delta }}{{\rm{\Omega }}}_{l}\sum _{m}^{3}\tfrac{{I}_{m}}{3}{n}^{(i)}{n}_{(j)}}{{\sum }_{l}{\rm{\Delta }}{{\rm{\Omega }}}_{l}{\sum }_{m}^{3}\tfrac{{I}_{m}}{3}},\end{eqnarray} \tag{ 48 }$$

where l indexes the triangulation faces, m indexes the vertices of each triangle, I_m is the frequency-integrated intensity at the mth vertex, and n⁽ⁱ⁾ is the unit vector for each vertex.

We now have ${R}^{0}{}_{\nu }$ in the coordinate frame and ${\pi }^{(i)}{}_{(j)}$ in the comoving frame; we want to recover a complete radiation stress tensor in the coordinate frame. The coordinate-frame four-momentum is related to the comoving frame stress–energy tensor by transformation matrices

$$\begin{eqnarray}&&{R}^{0}{}_{\nu }={e}_{(\mu )}{}^{0}{e}^{(\nu )}{}_{\nu }{R}^{(\mu )}{}_{(\nu )}.\end{eqnarray} \tag{ 49 }$$

Here, ${R}^{(i)}{}_{(j)}=-{R}^{(0)}{}_{(0)}{\pi }^{(i)}{}_{(j)}$ and ${R}^{(\nu )}{}_{(i)}=-{R}^{(i)}{}_{(\nu )}$, so we can expand the right-hand side to recover a system of four linear equations for the unknowns ${R}^{(0)}{}_{(\nu )}$. By inverting a 4 × 4 matrix, we recover ${R}^{(0)}{}_{(\nu )}$, which in turn yields ${R}^{(\mu )}{}_{(\nu )}$ through ${\pi }^{(i)}{}_{(j)}$, which can then be transformed back to ${R}^{\mu }{}_{\nu }$.

An initial set of samples will stream out of outflow boundary conditions on the light-crossing time of the simulation volume. Evidently, samples must be replenished in such situations. In addition, the number of samples per zone will fluctuate in multizone problems, especially when grids are irregular or gradients in Lorentz factor are appreciable, and we wish to ensure that the number of samples per zone never becomes too high or too low. Such a procedure also allows for dynamically controlling sample resolution: sensitive regions of integrations may require more samples than average.

Our approach to sample refinement and derefinement is simple and not unique. We first impose a desired number of samples per zone, which can in general depend on the local MHD or radiation properties. We then use this number of samples to define maximum and minimum solid angles for triangles on the unit sphere.

If our triangulation results in a triangle with a solid angle above the maximum, we create a sample nearly at the center of this face and randomly positioned inside the zone with a specific intensity that is the average of the triangle's vertices. Samples are not created exactly at the face center, to avoid providing the convex hull routine with two colocated vertices.

If our triangulation results in a triangle with a solid angle below the minimum, we identify the shortest edge of that triangle and randomly mark one of the two samples that form that edge for deletion.

Radiation interacts with plasma via emission, absorption, and scattering. For samples, these interactions are processed deterministically along characteristics. For the moment sector, we use sample intensities to frequency-average the opacities. Scattering is evaluated by using samples to construct an angle-averaged comoving spectrum, and evolving this spectrum with a scattering kernel. The change in four-momentum in the radiation moments determines the change in four-momentum of the fluid.

At large optical depths and/or small ratios of gas-to-radiation pressure, the characteristic timescale for four-momentum exchange between fluid and radiation can become much shorter than the global simulation time step, and the problem becomes stiff. Strong coupling can be stabilized with implicit methods. For time-dependent radiation transport, implicit methods can be applied on a per-zone or semi-implicit (i.e., the scheme remains explicit in spatial fluxes) basis, a major computational efficiency over global implicit solves. However, for intensities discretized in frequency and solid angle, such a semi-implicit solve could still require inverting a large matrix.

Instead, we adopt the "inner" and "outer" loop approach of Skinner et al. (2019), in which one rootfinds over only one or a few nonlinear equations in an outer loop (indexed by k), and in each step of that rootfind updates intensities in a semi-implicit fashion using the most recent (kth) value for the gas temperature to evaluate emissivities and absorptivities. We initialize this procedure by angle-averaging our sample intensities at time step n, where we have used linear interpolation to shift the comoving sample intensities ${I}_{\nu }^{n}$ onto a common frequency grid in the comoving frame:

$$\begin{eqnarray}&&{{ \mathcal I }}_{\nu }^{n}=\displaystyle \frac{{\sum }_{l}{\left({\rm{\Delta }}{\rm{\Omega }}\right)}_{l}{I}_{\nu }^{n}}{{\sum }_{l}{\left({\rm{\Delta }}{\rm{\Omega }}\right)}_{l}}.\end{eqnarray} \tag{ 50 }$$

3.8.1. Inner Loop

Our inner loop is composed of two steps. First, we use the Kompaneets equation to evaluate ${{ \mathcal J }}_{\nu }^{{\rm{s}},\,k}$ (Equation (32)) and ${\rm{\Delta }}{u}_{{\rm{r}}}^{{\rm{s}},\,k}$, the change in the comoving radiation energy density due to inelastic scattering, calculated as

$$\begin{eqnarray}&&{\rm{\Delta }}{u}_{{\rm{r}}}^{{\rm{s}},\,k}=\displaystyle \frac{4\pi }{c}{ \mathcal C }\sum _{i}\left({n}_{i}^{n+1}-{n}_{i}^{n}\right){\nu }_{i}{\rm{\Delta }}\mathrm{log}\nu .\end{eqnarray} \tag{ 51 }$$

In general, the integral of the ${{ \mathcal I }}_{\nu }$ evaluated from the samples will not correspond to the moment sector's ${u}_{{\rm{r}}}$. Rather than normalizing ${{ \mathcal I }}_{\nu }^{n}$ to recover ${u}_{{\rm{r}}}^{n}$, we normalize ${\rm{\Delta }}{u}_{{\rm{r}}}^{{\rm{s}}}$ by the ratio of the energy density evaluated in the sample sector to that of the moment sector, avoiding issues with maintaining thermal spectra in the samples.

We adopt the numerical method of Chang & Cooper (1970) for solving the Kompaneets equation. Here, we briefly review this method. We discretize the equation for the evolution of photon occupation number n over $x=h\nu /{k}_{{\rm{B}}}{T}_{{\rm{e}}}$ as:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{n}_{i}^{n+1}-{n}_{i}^{n}}{{\rm{\Delta }}\tau } & = & {n}_{{\rm{e}}}{\sigma }_{{\rm{T}}}c{{\rm{\Theta }}}_{{\rm{e}}}\displaystyle \frac{1}{{x}_{i}^{3}}\displaystyle \frac{1}{{\rm{\Delta }}\mathrm{log}{x}_{i}}\\ & & \times \left({x}_{i+1/2}^{4}{F}_{j+1/2}^{* }-{x}_{i-1/2}^{4}{F}_{j-1/2}^{* }\right),\end{array}\end{eqnarray} \tag{ 52 }$$

where Δτ is the elapsed proper time of the fluid from n to n + 1 and

$$\begin{eqnarray}&&{F}_{i+1/2}^{* }=\displaystyle \frac{{n}_{i+1}^{n+1}-{n}_{i}^{n+1}}{{x}_{i+1}-{x}_{i}}+\left(1+{n}_{i+1/2}^{n}\right){n}_{i+1/2}^{n+1},\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{n}_{i+1/2}^{n}=(1-{\delta }_{i}){n}_{i+1}^{n}+{\delta }_{i}{n}_{i}^{n},\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{n}_{i+1/2}^{n+1}=(1-{\delta }_{i}){n}_{i+1}^{n+1}+{\delta }_{i}{n}_{i}^{n+1},\end{eqnarray} \tag{ 55 }$$

where the $0\leqslant {\delta }_{i}\leqslant 0.5$ are chosen to ensure stability using a quasi-equilibrium distribution n_eq,

$$\begin{eqnarray}&&{w}_{i}=\left(1+{n}_{\mathrm{eq},j}/2+{n}_{\mathrm{eq},j+1}/2\right)\left({x}_{i+1}-{x}_{i}\right),\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\delta }_{i}=\displaystyle \frac{1}{{w}_{i}}-\displaystyle \frac{1}{\exp \left({w}_{i}\right)-1}.\end{eqnarray} \tag{ 57 }$$

Note that ${F}_{i+1/2}^{* }$ contains both ${n}_{i}^{n+1}$ and ${n}_{i}^{n};$ this discretization is not fully implicit, a potential source of instability—especially where n is much greater than unity (for example, scattering of lines with brightness temperatures much greater than the fluid temperature). We enforce ${F}^{* }=0$ at the boundaries, thereby conserving photon number to machine precision. After each update, we enforce a floor such that ${n}_{i}^{n+1}\geqslant {10}^{-100}$.

We solve this tridiagonal system of equations in ${ \mathcal O }({N}_{\mathrm{bin}})$ time, where N_bin is the number of frequency bins. This semi-implicit method can occasionally fail. We identify such situations by measuring the change in photon number density; if this quantity varies fractionally over a step by more than 10⁻⁸, we repeat the calculation without the nonlinear terms, yielding a fully implicit solution. One could also use an iterative method to solve the original equation in a fully implicit way; evidently, the nonlinear terms are important when the semi-implicit method fails. However, we adopt our simpler approach to ensure stability and avoid nonlinear multidimensional rootfinding, at the potential cost of approaching Wien rather than Bose–Einstein photon distributions in scattering-dominated media. Other methods for discretizing the Kompaneets equation (e.g., Larsen et al. 1985) may be less susceptible to this issue.

Another numerical difficulty is that calculation of the δ_i requires evaluation of a quasi-equilibrium solution. Here, we follow Chang & Cooper (1970) and choose a Bose–Einstein distribution:

$$\begin{eqnarray}&&{n}_{\mathrm{eq}}=\displaystyle \frac{1}{{{Ce}}^{x}-1},\end{eqnarray} \tag{ 58 }$$

where C, related to the photon chemical potential, is evaluated such that n_eq and n yield the same photon number density. Here, n_eq is singular at a point x > 0 for C < 1, which occurs when the photon number is greater than that of a blackbody at temperature T_e. We avoid this difficulty by enforcing C ≥ 2 (i.e., n_eq ≤ 1 everywhere) when calculating the δ_i. We have not found this trick to damage stability or accuracy.

In the second step, we compute angle- and frequency-averaged opacities and integrated emissivities. In the kth iteration of the outer loop, we approximate the updated angle-averaged spectrum via a backward Euler discretization:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{ \mathcal I }}_{\nu }}{{\nu }^{3}}\right)}^{k}=\displaystyle \frac{\tfrac{{{ \mathcal J }}_{\nu }^{{\rm{a}},\,k}+{{ \mathcal J }}_{\nu }^{{\rm{s}},\,k}}{{\nu }^{2}}+\tfrac{{\left({{ \mathcal I }}_{\nu }/{\nu }^{3}\right)}^{n}}{{\rm{\Delta }}s}}{\nu ({\widetilde{\alpha }}_{\nu }^{{\rm{a}},\,k}+{\widetilde{\alpha }}_{\nu }^{{\rm{s}},\,k})+\tfrac{1}{{\rm{\Delta }}s}},\end{eqnarray} \tag{ 59 }$$

where ${\widetilde{\alpha }}_{\nu }$ is the intensity-weighted, angle-averaged opacity at frequency ν. With these updated intensities in hand, we compute the frequency- and angle-averaged opacities:

$$\begin{eqnarray}&&{\overline{\alpha }}^{{\rm{a}},k}=\displaystyle \frac{\int d{\rm{\Omega }}\sum \nu {\rm{\Delta }}\mathrm{log}\nu {I}_{\nu }^{k}{\alpha }_{\nu }^{{\rm{a}},\,k}}{\int d{\rm{\Omega }}\sum \nu {\rm{\Delta }}\mathrm{log}\nu {I}_{\nu }^{k}},\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{\overline{\alpha }}^{{\rm{s}},k}=\displaystyle \frac{\int d{\rm{\Omega }}\sum \nu {\rm{\Delta }}\mathrm{log}\nu {I}_{\nu }^{k}{\alpha }_{\nu }^{{\rm{s}},\,k}}{\int d{\rm{\Omega }}\sum \nu {\rm{\Delta }}\mathrm{log}\nu {I}_{\nu }^{k}},\end{eqnarray} \tag{ 61 }$$

where the summation is over frequency bins. We also integrate the emissivity in a similar fashion, rather than analytically, so that truncation error in only absorption coefficient integration does not prevent the fluid from thermalizing effectively for finite numbers of frequency bins:

$$\begin{eqnarray}&&{J}^{{\rm{a}},k}=\int d{\rm{\Omega }}\sum \nu {\rm{\Delta }}\mathrm{log}\nu {j}_{\nu }^{{\rm{a}},\,k}.\end{eqnarray} \tag{ 62 }$$

For clarity, we left these expressions in terms of integrals over solid angle. In practice, we employ the same quadrature rule we described previously.

3.8.2. Outer Loop

In GRRMHD, the four-force update is four coupled equations. Previous authors have evaluated this in a fully nonlinear way with a 4D rootfind (Roedig et al. 2012; Sa̧dowski et al. 2013; McKinney et al. 2014) or a linearized 4D solve (Foucart et al. 2015). Here, we adopt a more efficient approach, performing a 1D nonlinear solve for the gas energy density (e.g., Skinner et al. 2019) in the comoving frame. We then construct an entire four-force G_(μ).

For our 1D rootfind, we write the gas energy update implicitly, and solve iteratively for ${u}_{{\rm{g}}}^{n+1}$ with a secant method:

$$\begin{eqnarray}&&{u}_{{\rm{g}}}^{n+1}-{u}_{{\rm{g}}}^{n}=-{u}_{{\rm{r}}}^{n+1}\left({\left(\displaystyle \frac{{{ \mathcal I }}_{\nu }}{{\nu }^{3}}\right)}^{n+1},{u}_{{\rm{g}}}^{n+1}\right)+{u}_{{\rm{r}}}^{n},\end{eqnarray} \tag{ 63 }$$

where ${u}_{{\rm{r}}}={R}^{(0)}{}_{(0)}$. In the ${k}^{\mathrm{th}}$ iteration, we compute ${u}_{{\rm{r}}}^{k+1}$, our newest estimate of ${u}_{{\rm{r}}}^{n+1}$, from Equation (29) as

$$\begin{eqnarray}&&{u}_{{\rm{r}}}^{k+1}=\displaystyle \frac{{u}_{{\rm{r}}}^{n}+{\rm{\Delta }}\tau {J}^{{\rm{a}},k}}{1+{\rm{\Delta }}\tau {\alpha }^{{\rm{a}},k}}+{\rm{\Delta }}{u}_{{\rm{r}}}^{{\rm{s}},\,k}.\end{eqnarray} \tag{ 64 }$$

Next, we use Equation (63) to evaluate a residual as

$$\begin{eqnarray}&&r=({u}_{{\rm{g}}}^{k}-{u}_{{\rm{g}}}^{n})+({u}_{{\rm{r}}}^{k+1}-{u}_{{\rm{r}}}^{n}),\end{eqnarray} \tag{ 65 }$$

and use the secant formula to compute our next guess for the gas energy, ${u}_{{\rm{g}}}^{k+1}$. Convergence is reached when a sufficiently small residual is computed, at which point we set ${u}_{{\rm{g}}}^{n+1}={u}_{{\rm{g}}}^{k+1}$. This implicit approach maintains stability when energy exchange is significant over a time step ${\rm{\Delta }}\tau$. The secant method occasionally fails. In such cases, we repeat the rootfind using bisection.

Once ${u}_{{\rm{g}}}^{n+1}$ is found, the radiation four-force is then calculated as

$$\begin{eqnarray}&&{G}_{(0)}^{n+1}=-\displaystyle \frac{{u}_{{\rm{g}}}^{n+1}-{u}_{{\rm{g}}}^{n}}{{\rm{\Delta }}\tau },\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&{G}_{(i)}^{n+1}=\displaystyle \frac{\left({\left({R}^{(0)}{}_{(i)}\right)}^{n+1}-{\left({R}^{(0)}{}_{(i)}\right)}^{n}\right)}{{\rm{\Delta }}\tau },\end{eqnarray} \tag{ 67 }$$

where ${R}^{(0)}{}_{(i)}$ is updated via backward Euler as

$$\begin{eqnarray}&&{\left({R}^{(0)}{}_{(i)}\right)}^{n+1}\equiv \displaystyle \frac{{\left({R}^{(0)}{}_{(i)}\right)}^{n}}{1+{\rm{\Delta }}\tau \left({\overline{\alpha }}^{{\rm{a}},n+1}+{\overline{\alpha }}^{{\rm{s}},n+1}\right)}.\end{eqnarray} \tag{ 68 }$$

The four-force is then transformed to the lab frame and applied to ${T}^{0}{}_{\nu }$ and ${R}^{0}{}_{\nu }$ using Equations (20) and (21). We also update the individual sample intensities in a backward Euler manner:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right)}^{n+1}=\displaystyle \frac{{j}_{\nu }^{{\rm{a}},n+1}/{\nu }^{2}+{{ \mathcal J }}_{\nu }^{{\rm{s}},n+1}/{\nu }^{2}+{\left({I}_{\nu }/{\nu }^{3}\right)}^{n}/{\rm{\Delta }}s}{\nu {\alpha }_{\nu }^{{\rm{a}},n+1}+\nu {\alpha }_{\nu }^{{\rm{s}},n+1}+1/{\rm{\Delta }}s}.\end{eqnarray} \tag{ 69 }$$

From our experiments, this approach shows excellent stability. Essentially, we exploit the different forms of G₍₀₎ and G_(i) in the comoving frame using an operator split; while G₍₀₎ is the difference of an emission coefficient and an absorption coefficient, both of which can be large, three-momentum exchange has only one term: ${G}_{(i)}=-\overline{\alpha }{R}^{(0)}{}_{(i)}$. Thus, while we rely on a numerical rootfind to evaluate the energy exchange, we can analytically evaluate the implicit momentum update. Since these interaction terms are typically treated with backward Euler for stability, this simple operator split does not degrade the temporal order of accuracy of the scheme.

Our method for time integration largely parallels the second-order midpoint scheme in harm (Gammie et al. 2003). While we could process the radiation subsystem separately in a first-order fashion (e.g., Jiang et al. 2014a; Ryan et al. 2015; Foucart 2018), i.e., one source term evaluation per time step, we have found improved performance on transport tests when using a second-order-in-time update to the advective fluxes of the radiation stress–energy tensor. The midpoint method also allows larger CFL numbers than an Euler step, so it is unclear whether we actually increase the cost of our scheme by going to second-order accuracy in time for advective fluxes. However, source term updates are still first-order in time. We typically use CFL = 0.8−0.9, independent of optical depth. Regardless of order of temporal accuracy, it is crucial, as we do here, to process in advance radiation interactions for any radiation moments used to source advective fluxes, in order to correctly recover diffusion speeds in optically thick media.

Here, we enumerate a complete time step from n to n + 1. Tildes denote which quantities have been updated by radiation interactions at their current n.

• 1.  
  Advect MHD conserved variables, ${\widetilde{U}}_{{\rm{g}}}^{n}\to {U}_{{\rm{g}}}^{n+1/2}$ sourced by ${\widetilde{U}}_{{\rm{g}}}^{n}$ (Section 3.1), and apply boundary conditions.
• 2.  
  Advect radiation moments, ${\left({\widetilde{R}}^{0}{}_{\nu }\right)}^{n}\to {\left({R}^{0}{}_{\nu }\right)}^{n+1/2}$ sourced by ${\left({\widetilde{R}}^{\mu }{}_{\nu }\right)}^{n}$ (Section 3.2), and apply boundary conditions.
• 3.  
  Calculate angle-averaged comoving intensity ${{ \mathcal I }}_{\nu }^{n}$ (Equation (50)), to be used for scattering in both Steps 7 and 14.
• 4.  
  Push samples along geodesics, ${({x}^{i},{n}_{\mu })}^{n}\to {({x}^{i},{n}_{\mu })}^{n+1/2}$ and ${\left({\widetilde{I}}_{\nu }/{\nu }^{3}\right)}^{n}\to {\left({I}_{\nu }/{\nu }^{3}\right)}^{n+1/2}$ (Section 3.3), and apply boundary conditions.
• 5.  
  Boost samples to the fluid frame (Section 3.4) and construct triangulations in each zone (Section 3.5).
• 6.  
  Integrate sample intensities over frequency and solid angle using triangles to evaluate ${\left({\pi }^{(i)}{}_{(j)}\right)}^{n+1/2}$ and then ${\left({R}^{i}{}_{j}\right)}^{n+1/2}$ (Section 3.6).
• 7.  
  Calculate radiation four-force at n+1/2, apply to conserved MHD and radiation quantities, ${U}_{{\rm{g}}}^{n+1/2}\to {\widetilde{U}}_{{\rm{g}}}^{n+1/2}$ and ${\left({R}^{0}{}_{\mu }\right)}^{n+1/2}\to {\left({\widetilde{R}}^{0}{}_{\mu }\right)}^{n+1/2}$, and sample intensities ${\left({I}_{\nu }/{\nu }^{3}\right)}^{n+1/2}\to {\left({\widetilde{I}}_{\nu }/{\nu }^{3}\right)}^{n+1/2}$ (Section 3.8).
• 8.  
  Recalculate ${\left({\widetilde{\pi }}^{(i)}{}_{(j)}\right)}^{n+1/2}$ with the ${\widetilde{I}}_{\nu }^{n+1/2}$ using the same triangulation, and use to evaluate ${\left({\widetilde{R}}^{i}{}_{j}\right)}^{n+1/2}$
• 9.  
  Advect MHD conserved variables, ${\widetilde{U}}_{{\rm{g}}}^{n}\to {U}_{{\rm{g}}}^{n+1}$ sourced by ${\widetilde{U}}_{{\rm{g}}}^{n+1/2}$ (Section 3.1), and apply sources to sample boundary conditions.
• 10.  
  Advect radiation moments, ${\left({\widetilde{R}}^{0}{}_{\nu }\right)}^{n}\to {\left({R}^{0}{}_{\nu }\right)}^{n+1}$ sourced by ${\left({\widetilde{R}}^{\mu }{}_{\nu }\right)}^{n+1/2}$ (Section 3.2), and apply boundary conditions.
• 11.  
  Push samples along geodesics, ${({x}^{i},{n}_{\mu })}^{n+1/2}\,\to {({x}^{i},{n}_{\mu })}^{n+1}$ and ${\left({\widetilde{I}}_{\nu }/{\nu }^{3}\right)}^{n+1/2}\to {\left({I}_{\nu }/{\nu }^{3}\right)}^{n}$ (Section 3.3), and apply boundary conditions.
• 12.  
  Boost samples to the fluid frame (Section 3.4) and construct triangulations in each zone (Section 3.5).
• 13.  
  Integrate sample intensities over frequency and solid angle using triangles to evaluate Eddington tensor ${\left({\pi }^{(i)}{}_{(j)}\right)}^{n+1}$ and ${\left({R}^{i}{}_{j}\right)}^{n+1}$ (Section 3.6).
• 14.  
  Calculate radiation four-force at n + 1, apply to conserved MHD and radiation quantities, ${U}_{{\rm{g}}}^{n+1}\to {\widetilde{U}}_{{\rm{g}}}^{n+1}$ and ${\left({R}^{0}{}_{\mu }\right)}^{n+1}\to {\left({\widetilde{R}}^{0}{}_{\mu }\right)}^{n+1}$, and apply sources to sample intensities ${\left({I}_{\nu }/{\nu }^{3}\right)}^{n+1}\to {\left({\widetilde{I}}_{\nu }/{\nu }^{3}\right)}^{n+1}$ (Section 3.8).
• 15.  
  Recalculate ${\left({\widetilde{\pi }}^{(i)}{}_{(j)}\right)}^{n+1}$ with the ${\widetilde{I}}_{\nu }^{n+1}$ using the same triangulation, and use to evaluate ${\left({\widetilde{R}}^{i}{}_{j}\right)}^{n+1}$

Consider a hohlraum boundary condition at x = 0 and temperature T such that I_ν = B_ν(T), and a vacuum for x > 0. The radiation energy density at the boundary is ${u}_{{\rm{r}},0}={a}_{{\rm{r}}}{T}^{4}$. Radiation will propagate in the positive x direction.

The time-dependent analytic solution is evaluated by sending characteristics backward in time over all θ from position x ($\theta =0$ corresponds to the $+x$ direction) and determining whether they reach the hohlraum boundary by t = 0. The specific intensity is then

$$\begin{eqnarray}{I}_{\nu }=\left\{\begin{array}{ll}{B}_{\nu } & \theta \lt {\theta }_{\max }\\ 0 & \mathrm{else}\end{array}\right.,\end{eqnarray} \tag{ 74 }$$

where ${\theta }_{\max }={\cos }^{-1}\left(x/{ct}\right)$, the radiation energy density is

$$\begin{eqnarray}&&{u}_{{\rm{r}}}=2\pi {\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{{\theta }_{\max }}{B}_{\nu }\sin \theta d\theta d\nu ,\end{eqnarray} \tag{ 75 }$$

and the Eddington factor is

$$\begin{eqnarray}&&{f}^{(x)}{}_{(x)}=1/3\left[{\left(\displaystyle \frac{x}{{ct}}\right)}^{2}+\displaystyle \frac{x}{{ct}}+1\right].\end{eqnarray} \tag{ 76 }$$

At late time, ${\theta }_{\max }\to \pi /2;$ while the intensity is far from thermal, the Eddington factor ${f}^{(x)}{}_{(x)}=1/3$, exactly the value for I_ν = B_ν. Although simple, this problem is analogous to physical situations encountered in radiation transport, in which a thermal object radiates into a tenuous atmosphere or ambient medium in nearly plane-parallel symmetry.

The hohlraum boundary condition in our simulation is enforced by setting the radiation moments to the value for a blackbody at temperature T. For MOCMC, we randomly distribute samples in the boundary zone at each time step, with I_ν = B_ν for these samples. However, because we are studying convergence on this problem, we want to avoid the discontinuity in radiation energy density at x = 0. Instead, we consider the boundary condition slightly away from x = 0, where only characteristics moving in the +x direction are thermal. Hence, our boundary condition is

$$\begin{eqnarray}{I}_{\nu }=\left\{\begin{array}{ll}{B}_{\nu } & \theta \leqslant \pi /2\\ 0 & \mathrm{else}\end{array}\right.\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}{R}^{\mu }{}_{\nu }=\left[\begin{array}{cccc}-{u}_{{\rm{r}},0}/2 & {u}_{{\rm{r}},0}/4 & 0 & 0\\ -{u}_{{\rm{r}},0}/4 & {u}_{{\rm{r}},0}/6 & 0 & 0\\ 0 & 0 & {u}_{{\rm{r}},0}/6 & 0\\ 0 & 0 & 0 & {u}_{{\rm{r}},0}/6\end{array}\right]\end{eqnarray} \tag{ 78 }$$

at x = 0. Note that Eddington and M1 produce similar results for this and an isotropic thermal boundary with I_ν = B_ν and the corresponding ${R}^{\mu }{}_{\nu }$. The right boundary condition is placed at large enough x to be causally disconnected from the simulation region shown.

Figure 2 compares Eddington, M1, and MOCMC closures on this test at t = 0.75 and t = 5. The assumption in M1, that there is a frame in which the radiation is isotropic, is not satisfied here. Additionally, ${f}^{(x)}{}_{(x)}=1/3$ at late time while f^(x) > 0, which is inconsistent with Equations (71) and (72). Even in a time-independent sense, M1 leads to an order unity error in the radiation energy density, unlike the Eddington closure.⁶ MOCMC accurately matches the analytic solution for both u_r and ${f}^{(x)}{}_{(x)}$.

Zoom In Zoom Out Reset image size

Convergence is shown in Figure 3. Evidently, we recover approximately first-order convergence in N_samp, unlike the ${N}_{\mathrm{samp}}^{-1/2}$ for a Monte Carlo method. This is a result of using the Delaunay triangulation to calculate ${\pi }^{(i)}{}_{(j)};$ essentially, the triangulation provides a more accurate estimate of the solid angle owned by each sample. We do not expect convergence better than first order on this test, which contains a discontinuity in momentum space.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now consider a 2D box with hohlraum boundary conditions on both the x = 0 and y = 0 boundaries.

We construct a solution for comparison via a simple Monte Carlo method. On a grid of positions (x, y) and time t, we sample characteristics uniformly on the sphere, propagate them back to t = 0, and ask whether they intersect a hohlraum boundary. Summing intensities over a solid angle in each zone gives u_r.

Unlike in the 1D hohlraum test, the modified boundaries with nonzero fluxes are no longer valid, so we set the boundary conditions to simply be thermal at x = 0 and y = 0. The results for Eddington, M1, and MOCMC compared to the semianalytic solution are shown in Figure 4.

This test shows that multidimensional effects are important to monitor when discriminating between transport algorithms. In particular, both moment methods exhibit significant interactions between the wavefronts from the two boundary conditions. These interactions lead to much larger radiation energy densities in the moment methods than are encountered in either MOCMC or the semianalytic solution. The MOCMC solution at 64 samples per zone also exhibits some radiation self-interaction due to truncation error in the Eddington tensor evaluation, but this self-interaction decreases with the number of samples and is already much lower than that seen in the moment methods.

Figure 5 shows the change in the MOCMC solution with increasing sample resolution from 16 to 128 samples per zone. While the boundary prescription we adopt in this multidimensional test makes formal convergence testing difficult, the MOCMC solution appears to approach the semianalytic solution in most of the domain as the sample resolution is increased.

Zoom In Zoom Out Reset image size

Thermalization provides a useful test of the code's ability to recover a basic feature of radiation hydrodynamics—equilibration between material and radiation temperatures. We repeat the bremsstrahlung thermalization test from Ryan et al. (2015) on a 3 × 3 grid of spatial zones, with a much larger initial gas temperature ${T}_{{\rm{g}},0}={10}^{9}\,{\rm{K}}$ and electron number density ${n}_{{\rm{e}}}=6\times {10}^{16}\,{\mathrm{cm}}^{-3}$. There is no radiation initially. The gas and radiation are allowed to proceed toward equilibrium. We compare to a frequency-dependent semianalytic solution in Figure 6.

Zoom In Zoom Out Reset image size

We also compare to frequency-integrated source terms i.e., we use a Planck mean opacity rather than an opacity averaged over the samples. Because the exponential cutoff in the emissivity shifts dramatically downward in frequency, the timescale to thermalize radiation that is emitted initially is very long, resulting in an undershoot of the gas temperature. This is not captured by a gray method, which cannot know the frequency distribution of radiation. Figure 7 shows the frequency distribution of radiation at different times during this test.

Zoom In Zoom Out Reset image size

The time evolution of a temperature perturbation in an otherwise uniform medium in radiative equilibrium with gray absorption coefficient α permits an exact solution under the assumption that the time-dependent terms in the radiative transfer equation are small (Spiegel 1957; Mihalas & Mihalas 1984). Note that this test assumes that light is transported much faster than the mode decay time; attempting to recreate this test in a slow-light code like ours will generally introduce some numerical diffusion through the Riemann solver in the radiation sector.

Unlike Ryan et al. (2015), here we hold α₀, T₀, and ${{ct}}_{\mathrm{RR}}/L$ constant. This problem requires that the relaxation rate be slow compared to the light-crossing time for the perturbation wavelength. We use 32 grid zones to simulate one wavelength.

We measure the relaxation rate by fitting the form of the linear solution to the numerical result after one e-folding time. We compare to the gas temperature only, because the perturbed radiation energy density is not trivial when the optical depth is finite. We show the performance across optical depth by comparing measured dispersion rates for Eddington, M1, and MOCMC to the analytic solutions, both for transport and for the Eddington closure, in Figure 8.

Zoom In Zoom Out Reset image size

We also show the analytic solution for the Eddington approximation. For the moderately optically thick limit, these closures produce similar errors, ≲20%. M1 does not distinguish between isotropic and anisotropic components, and the isotropic component dominates the four-momentum used to compute ${f}^{(x)}{}_{(x)}$.

The Eddington factor for the perturbation is

$$\begin{eqnarray}&&{f}^{(x)}{}_{(x)}\left(\tau \right)=\displaystyle \frac{\tau }{2\pi }\displaystyle \frac{1-\tfrac{\tau }{2\pi }{\tan }^{-1}\left(\tfrac{2\pi }{\tau }\right)}{{\tan }^{-1}\left(\tfrac{2\pi }{\tau }\right)};\end{eqnarray} \tag{ 79 }$$

for $\tau \to \infty$, ${f}^{(x)}{}_{(x)}\to 1/3;$ as expected, the Eddington approximation is appropriate. For $\tau \to 0$, however, ${f}^{(x)}{}_{(x)}\to 0$, whereas Eddington and M1 closure both assume ${f}^{(x)}{}_{(x)}\geqslant 1/3;$ the perturbation in the radiation field becomes oblate along the x axis. Eddington factors less than 1/3 also appear in radiative shocks (Jiang et al. 2014a). Clearly, standard moment closure models are invalid here, and yet such closures recover the correct dispersion relation at small τ. Mihalas & Mihalas (1984) point out that this is because, for $\tau \to 0$, relaxation is determined entirely by emission.

We repeat the setup from Ryan et al. (2015) of thermalization of soft photons due to Compton upscattering in a one-zone box. For this test, we use 32 samples per zone, although there is no angular structure, so there is no truncation error due to solid angle discretization (note the absence of noise in the solution). We use 100 frequency bins, logarithmically spaced from ${10}^{8}$ to ${10}^{20}\,\mathrm{Hz}$. We initialize the gas with electron and proton number densities ${n}_{{\rm{e}}}=2.5\times {10}^{17}\,{\mathrm{cm}}^{-3}$ and temperature ${T}_{{\rm{g}},0}=5\times {10}^{7}\,{\rm{K}}$. The radiation is initially monochromatic at frequency ${\nu }_{0}=3\times {10}^{16}\,\mathrm{Hz}$ and photon number density ${n}_{\gamma }=2.38\times {10}^{18}\,{\mathrm{cm}}^{-3}$.

The characteristic timescales in this problem are the mean time between scattering events, ${t}_{{\rm{s}}}=1/({n}_{{\rm{e}}}{\sigma }_{{\rm{T}}}c)\approx 2\,\times \,{10}^{-4}\,{\rm{s}}$, and the Comptonization time, i.e., the time between scatterings divided by the fractional energy transfer per scattering event, ${t}_{{\rm{C}}}={t}_{{\rm{s}}}{m}_{{\rm{e}}}{c}^{2}/({k}_{{\rm{B}}}{T}_{{\rm{e}}})\approx 2.4\,\times \,{10}^{-2}\,{\rm{s}}$.

We calculate the equilibrium temperature with conservation of energy and photon number. Unlike Ryan et al. (2015), here we assume that the final photon distribution is Bose–Einstein rather than Wien; the final temperature T_f is found by solving

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{2{n}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{g}},0}}{\gamma -1}+h{\nu }_{0}{n}_{\gamma } & = & \displaystyle \frac{2{n}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{f}}}}{\gamma -1}\\ & & +\,\displaystyle \frac{48\pi {k}_{{\rm{B}}}^{4}{T}_{{\rm{f}}}^{4}}{{c}^{3}{h}^{3}}{\mathrm{Li}}_{4}\left(\exp \left(\displaystyle \frac{-\mu }{{k}_{{\rm{B}}}{T}_{{\rm{f}}}}\right)\right),\end{array}\end{eqnarray} \tag{ 80 }$$

where μ is the chemical potential of the photons and ${\mathrm{Li}}_{s}\left(z\right)$ is the polylogarithm. For this test, the photon occupation number remains much less than one, and T_f closely agrees with the value calculated assuming a final Wien distribution for the photons.

We perform this test both with an initially isotropic distribution of radiation, and with an anisotropic distribution in photon number, ${n}_{\gamma ,\mathrm{aniso}}={n}_{\gamma }(1+{n}^{(x)}),$ with 64 samples. The gas and radiation temperature (the radiation temperature is evaluated assuming a Boltzmann distribution) evolutions for both cases are shown in Figure 9. The gas and radiation approach thermal equilibrium on a timescale $\approx {t}_{{\rm{C}}}$, as expected. A similar result was obtained for the same parameters with the bhlight code (Ryan et al. 2015). The evolution of the intensity in angle is shown in Figure 10. Because the fractional energy change of each scattering event is small, the radiation isotropizes rapidly compared to the timescale for thermalization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To test the performance of our scattering treatment, as well as the behavior of MOCMC in the diffusion regime, we consider diffusion of a Gaussian pulse in a static medium optically thick to Thomson scattering in 1D on a domain $x\in \left[-L/2,L/2\right]$. Starting at t₀, the analytic solution for the radiation energy density (with a 0.01% background) in the diffusion regime is

$$\begin{eqnarray}&&{u}_{{\rm{r}}}={u}_{{\rm{r}},0}\left({10}^{-4}+\sqrt{\displaystyle \frac{{t}_{0}}{t}}\exp \left(\displaystyle \frac{-{x}^{2}}{4{Dct}}\right)\right),\end{eqnarray} \tag{ 81 }$$

where the diffusion coefficient $D=1/(3{\sigma }_{{\rm{T}}}{n}_{{\rm{e}}})$ and ${u}_{{\rm{r}},0}$ is the maximum radiation energy density at the initial time t = t₀. We set n_e such that the optical depth over the domain is τ = 10⁴. We evolve this system to t = 50L/c; the MOCMC solution is shown in Figure 11. The solution shows good agreement; in particular, the pulse diffuses much more slowly than the rate at which samples traverse zones.

Zoom In Zoom Out Reset image size

To study dynamic diffusion, we add a lab-frame velocity to the fluid such that the radiation is advected with the flow more rapidly than it diffuses. This is a powerful test of accuracy for a radiation hydrodynamics method. In particular, this test can be a challenge for ${ \mathcal O }\left(v/c\right)$ methods when care is not taken in truncating the equations in a way appropriate for all regimes (Krumholz et al. 2007). This is captured naturally by covariant methods like MOCMC.

We adopt the domain and initial conditions from the previous test, except that the pulse is now initially centered at x/L = −0.25. We set the fluid speed to β ≡ v/c = 0.1, and we evolve the system to t_f = 5L/c. As in the previous test, we set τ = 10⁴; at t_f, there is some diffusion of the pulse. Figure 12 shows the MOCMC solution, along with initial and final analytic solutions.

Zoom In Zoom Out Reset image size

We consider a 1D box of length L with gas and radiation initially in thermal equilibrium at temperature T₀ = 10⁷ K. The gas and radiation interact through a gray absorption opacity. As this system evolves, noise in the Eddington tensor will lead to noise in the radiation four-momentum, which in turn will couple to the gas energy density. We use this test to measure the noise in the fluid when varying the optical depth across the box τ and the gas-to-radiation pressure ratio β_r. Note that noise-free methods like analytic moment closures satisfy this test trivially.

We parameterize the stability by $\langle | {T}_{{\rm{g}}}-{T}_{0}| \rangle /{T}_{0}$, where T_g is the gas temperature and $\langle \cdot \rangle$ is the average over the domain. We set β_r and τ by setting the gas density and the gray opacity κ. We run each realization for t = 10L/c, which appears to lead to saturated gas temperature errors. We consider a range of optical depths and gas-to-radiation pressures, $\tau \in [{10}^{-2},{10}^{2}]$, and ${\beta }_{{\rm{r}}}\in [{10}^{-5},{10}^{2}]$. The mean error at final time for each of these realizations is shown in Figure 13. Evidently, our method is stable for every combination of parameters we consider, and mean gas temperature errors remain good $(\lesssim 0.5 \% )$ even in the most extreme cases.

Zoom In Zoom Out Reset image size

Nonetheless, noise in the gas temperature does increase with optical depth and radiation to gas pressure. This test already considers large optical depths and very large radiation to gas pressure ratios, but even more extreme parameter regimes may obtain in problems of astrophysical relevance. One option would be to interpolate toward the (noise-free) Eddington closure in very optically thick regions, although this could cause issues if there also exists some radiation at frequencies that are optically thin.

We revisit relativistic radiation MHD linear modes (Jiang et al. 2012; Sa̧dowski et al. 2014; Ryan et al. 2015). These are derived assuming the Eddington closure, and so we focus on the optically thick regime here. Solutions are generated with a symbolic linear modes package (Chandra et al. 2017). We specialize to optically thick radiation-modified fast magnetosonic modes at different gas-to-radiation pressures β_r. We do not refine or derefine samples for this test. We construct eigenmodes of the form $\delta \sim \exp \left(\omega t+{ikx}\right)$ for variation in ${\boldsymbol{P}}=({\rho }_{0}+\delta \rho$, ${u}_{{\rm{g}},0}+\delta {u}_{{\rm{g}}}$, $\delta {u}^{1}$, δu², ${B}_{0}^{1}$, ${B}_{0}^{2}+\delta {B}^{2}$, ${u}_{{\rm{r}},0}+\delta {u}_{{\rm{r}}}$, $\delta {F}^{1}$, $\delta {F}^{2})$. Optical depth per wavelength $\tau =20$, divided evenly between gray absorption opacity ${\kappa }^{{\rm{a}}}$ and gray scattering opacity ${\kappa }^{{\rm{s}}}$, for wavenumber $k=2\pi$. Modes are simulated with an amplitude $\delta ={10}^{-4}$. The background equilibrium is ${\rho }_{0}=1$, ${u}_{0}^{1}={u}_{0}^{2}={F}_{0}^{1}={F}_{0}^{2}=0$, ${E}_{0}={a}_{{\rm{r}}}{({P}_{0}/{\rho }_{0})}^{4}$, and ${u}_{{\rm{g}},0}$ and ${B}_{0}^{1}={B}_{0}^{2}$ are determined by plasma ${\beta }_{{\rm{m}}}=1$ and gas-to-radiation pressure ${\beta }_{{\rm{r}}}=(10,1,0.1)$. The adiabatic index $\gamma =5/3$. Units are such that ${k}_{{\rm{B}}}=c={a}_{{\rm{r}}}=1$. The three modes we consider are given in Table 1.

Each mode is simulated for one wave crossing time, $2\pi /| \mathrm{Imag}\left(\omega \right)|$. We study convergence. However, we have two resolution parameters: number of zones N¹, and number of samples per zone ${N}_{\mathrm{samp}}$, which introduce similar errors. We therefore approach convergence in two steps. First, we study convergence of the MOCMC code with the Eddington closure relative to the analytic solution with increasing N¹ (this subsystem has no ${N}_{\mathrm{samp}}$ dependence). Figure 14 shows the expected ${({N}^{1})}^{-1}$ convergence in all variables at large N¹ (these modes are optically thick, and our source term evaluations are first-order accurate in time. Next, we fix ${N}^{1}=256$, and reintroduce samples for computing the pressure tensor. We then study convergence of the full MOCMC solution with increasing ${N}_{\mathrm{samp}}$ relative to the analytic solution. When the truncation error is dominated by the integration of the Eddington tensor, MOCMC converges as $\sim {N}_{\mathrm{samp}}^{-2}$.

Zoom In Zoom Out Reset image size

Farris et al. (2008) introduced a method for calculating 1D relativistic radiation hydrodynamic waves in flat spacetime with an assumed Eddington closure and a gray absorption opacity, along with four example solutions that have frequently been reproduced with relativistic radiation hydrodynamics codes (Fragile et al. 2012; Roedig et al. 2012; Sa̧dowski et al. 2013; McKinney et al. 2014; Ryan et al. 2015). These four example solutions are: (Case 1) a gas pressure–dominated nonrelativistic strong shock, (Case 2) a gas pressure–dominated mildly relativistic strong shock, (Case 3) a gas pressure–dominated highly relativistic wave, and (Case 4) a radiation pressure–dominated, mildly relativistic wave. We adopt the parameters from Farris et al. (2008). Ryan et al. (2015) simulated Cases 1–3 with full transport. See also Ohsuga & Takahashi (2016) for another full transport method applied to this problem.

We initialize these problems as shocktubes and allow them to evolve to equilibrium inside the code. We perform these simulations with both the Eddington closure and the full MOCMC machinery. The Eddington closure provides the reference solution and tests part of our numerical framework against the analytic solution (Farris et al. 2008); the MOCMC solution, being full-transport, will not agree with the analytic solution on scales of an optical depth, which in all cases is approximately the scale of the interface structure.

For all cases, we use 800 spatial zones, and for the MOCMC solution, approximately 64 samples per zone after refinement. Each sample carries 50 frequency bins. The left and right initial interface states are enforced at the boundary in the fluid and radiation moment variables, while the pressure tensor is taken to be Eddington and the samples are thermal and uniformly distributed in solid angle in the comoving frame. The four waves are shown, respectively, in Figures 15–18. The particle noise in the MOCMC representation of the Eddington tensor is mostly small and does not prevent the code from being stable or accurate in any case, even when radiation pressure is dominant. The most significant pathology is when the Lorentz factor changes dramatically across the wave (Case 3); the resulting beaming of samples leads to poor sampling of solid angle in the comoving frame. The code then relies on the resampling procedure for resolving the pressure tensor.

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

We now consider radiation transport in curved spacetime. We essentially repeat the 2D flat space hohlraum test in the Kerr geometry for a $1\,{M}_{\odot }$ black hole with spin a = 0.9375, now with the radiating boundary condition a thin disk at the midplane from $6\,{GM}/{c}^{2}$ to $10\,{GM}/{c}^{2}$. The disk has the temperature profile of a thin disk (Novikov & Thorne 1973) with anomalous viscosity $\alpha =0.05$. The disk and atmosphere are static in the normal observer frame. This disk radiates into vacuum; at finite time, we measure the radiation energy density in the normal observer frame.

To construct a time-dependent solution to the equation of radiative transfer, we extend the procedure from 4.2 to a geometrically thin, optically thick disk in the Kerr geometry. The procedure is essentially unchanged—except now, instead of propagating rays backward in time along straight lines in flat space, we sample geodesics uniformly in the normal observer frame (e.g., Bardeen et al. 1972) and propagate them along geodesics until they either exit the outer radial boundary, cross the event horizon, or intersect with the thin disk. For geodesics that intersect the disk, we set their invariant intensity to the Planck function at the disk temperature, and integrate this over frequency in the originating normal observer frame to get the contribution to radiation energy density ${u}_{{\rm{r}},\mathrm{normal}}$ in that frame.

Figure 19 shows the results of this test, both for the full MOCMC code as well as Eddington and M1 closures, alongside the semianalytic solution. For the simulations, we adopt axisymmetry in modified Kerr–Schild (Gammie et al. 2003) coordinates with refinement parameter h = 0.3 and use 128 × 129 zones in X¹ and X², the odd number of zones in X² for symmetry about the midplane. The lower hemisphere is not shown. We set the outer radius at $20{GM}/{c}^{2}$. The thermal radiating disk boundary (we adopt a hemispheric approach similar to the 2D hohlraum test in Section 4.2) is enforced at the midplane every substep, both in the radiation moments and in the samples (i.e., the radiating disk is infinitely optically thick in our implementation).

Zoom In Zoom Out Reset image size

We repeat the isothermal pressure-supported atmosphere test in the Schwarzschild geometry close to the event horizon from Ryan et al. (2015). Reflecting spherical shells are placed at inner and outer radii ${r}_{\mathrm{in}}=3.5{GM}/{c}^{2}$ and ${r}_{\mathrm{out}}=20{GM}/{c}^{2}$, and gas between these shells is allowed to reach radiative equilibrium through a gray opacity κ. The radiation acts like a heat conduction, leading to a temperature profile that is isothermal modulo a redshift factor,

$$\begin{eqnarray}&&T\left(r\right)={T}_{\infty }/\sqrt{-{g}_{00}}\end{eqnarray} \tag{ 82 }$$

where T_∞ is the temperature at large radius. The solution is determined by this temperature profile and mechanical equilibrium, ${T}^{\mu r}{}_{;\mu }=0$. For a γ-law equation of state, the solution is evaluated by solving

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dR}}=-\displaystyle \frac{\left(\rho +\tfrac{\gamma }{\gamma -1}{P}_{{\rm{g}}}\right)}{{r}^{2}\left(1-\tfrac{2}{r}\right)}.\end{eqnarray} \tag{ 83 }$$

Because all rays propagate back to $t=-\infty$, ${I}_{\nu }/{\nu }^{3}={B}_{\nu }/{\nu }^{3}$ everywhere in the domain, and this is an exact solution of the equations of radiation hydrodynamics with full transport in curved spacetime. This problem can be cast in terms of three dimensionless parameters: (1) the ratio of the inner atmospheric scale height H to r_in, $H/{r}_{\mathrm{in}}={k}_{{\rm{B}}}{T}_{\mathrm{in}}{r}_{\mathrm{in}}/(\mu {m}_{{\rm{p}}}{GM})$ where μ is the mean molecular weight; (2) the ratio of gas-to-radiation pressure at the inner boundary, ${\beta }_{{\rm{r}}}=\mu {m}_{{\rm{p}}}{a}_{{\rm{r}}}{T}_{\mathrm{in}}^{3}/(3{\rho }_{\mathrm{in}}{k}_{{\rm{B}}})$; and (3) the optical depth across the domain $\tau =\kappa {\rho }_{\mathrm{in}}({r}_{\mathrm{out}}-{r}_{\mathrm{in}})$. We set the black hole mass M = M_⊙, $H/{r}_{\mathrm{in}}=1.60$, ${\beta }_{{\rm{r}}}=43.5$, and τ = 5. We set up the simulation in 1D, with 128 grid zones. The exact solution is enforced at the boundaries. We run for $t=500{GM}/{c}^{3}$. The result is shown in Figure 20.

Zoom In Zoom Out Reset image size