bhlight: GENERAL RELATIVISTIC RADIATION MAGNETOHYDRODYNAMICS WITH MONTE CARLO TRANSPORT

We assume particle number conservation, which in a coordinate basis is

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

where ρ₀ is the comoving frame rest mass density and u^μ is the fluid four-velocity.

Energy and momentum conservation for the coupled fluid and radiation system are given by

$$\begin{eqnarray}&&{\left({T}_{\nu }^{\mu }+{R}_{\nu }^{\mu }\right)}_{;\mu }=0,\end{eqnarray} \tag{ 2 }$$

where ${T}_{\nu }^{\mu }$ is the magnetohydrodynamic stress–energy tensor, and ${R}_{\nu }^{\mu }$ is the radiation stress–energy tensor (not to be confused with the Ricci tensor). In a coordinate basis, Equation (2) becomes

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

where the radiation four-force density

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

In the ideal MHD limit u_μ F^μν = 0 (F^μν ≡ electromagnetic field tensor), and one can show that

$$\begin{eqnarray}&&{T}_{\nu }^{\mu }=({\rho }_{0}+u+P+{b}^{2}){u}^{\mu }{u}_{\nu }+(P+\displaystyle \frac{1}{2}{b}^{2}){g}_{\nu }^{\mu }-{b}^{\mu }{b}_{\nu },\end{eqnarray} \tag{ 5 }$$

where P ≡ fluid pressure, u ≡ fluid internal energy density, and b² = b^μb_μ, with b^μ the magnetic field four-vector,

$$\begin{eqnarray}&&{b}^{\mu }\equiv \displaystyle \frac{1}{2}{\epsilon }^{\mu \nu \kappa \lambda }{u}_{\nu }{F}_{\lambda \kappa },\end{eqnarray} \tag{ 6 }$$

and ${\epsilon }^{\mu \nu \kappa \lambda }\equiv -[\mu \nu \kappa \lambda ]/\sqrt{-g}$ is the Levi–Civita tensor.

Evidently b^μu_μ = 0, so b^μ has only three degrees of freedom, expressed as Bⁱ ≡ ^*F^it, where

$$\begin{eqnarray}&&{}^{*}{F}^{\mu \nu }=\displaystyle \frac{1}{2}{\epsilon }_{\mu \nu \kappa \lambda }{F}^{\kappa \lambda }={b}^{\mu }{u}^{\nu }-{b}^{\nu }{u}^{\mu }.\end{eqnarray} \tag{ 7 }$$

Then

$$\begin{eqnarray}&&{b}^{t}={B}^{i}{u}^{\mu }{g}_{i\mu },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{b}^{i}=\displaystyle \frac{{B}^{i}+{b}^{t}{u}^{i}}{{u}^{t}}.\end{eqnarray} \tag{ 9 }$$

The magnetic field evolution is determined by

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

subject to the no-monopoles constraint

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

The equation of state is

$$\begin{eqnarray}&&P=(\gamma -1)u.\end{eqnarray} \tag{ 12 }$$

To summarize, the governing equations for the fluid evolution are Equations (1), (3), and (10), together with (11) and (12).

The radiation field consists of photons with wave four-vector k^μ and momentum ${p}^{\mu }={\rm{\hslash }}{k}^{\mu }$. The photons follow geodesics, with

$$\begin{eqnarray}&&\displaystyle \frac{{{dx}}^{\mu }}{d\lambda }={k}^{\mu }\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{{dk}}^{\mu }}{d\lambda }=-{{\rm{\Gamma }}}_{\mu \nu }^{\lambda }{k}^{\mu }{k}^{\nu },\end{eqnarray} \tag{ 14 }$$

where ${{\rm{\Gamma }}}_{\mu \nu }^{\lambda }$ is the connection and λ is an affine parameter along the geodesic. We assume that plasma dispersion effects are negligible, so photons travel on null geodesics, k_μ k^μ = 0. The frequency of a photon in a frame with four-velocity u^μ is ω = −k^μu_μ ($\nu \equiv \omega /(2\pi )$).

In nonrelativistic radiative transfer one describes the radiation field with the specific intensity I_ν (here and throughout we ignore polarization), which is frame-dependent. However, I_ν/ν³ ∝ f_R where f_R is the radiation distribution function

$$\begin{eqnarray}&&{f}_{{\rm{R}}}({x}^{\mu },{p}_{i})=\displaystyle \frac{{dN}}{{d}^{3}{{xd}}^{3}p},\end{eqnarray} \tag{ 15 }$$

where d³p = dp₁dp₂dp₃. Because dN, ${d}^{3}{{xp}}^{t}\sqrt{-g},$ and ${d}^{3}p/({p}^{t}\sqrt{-g})$ are invariant, f_R is also invariant.

The evolution of f_R is given by the Boltzmann equation

$$\begin{eqnarray}&&\displaystyle \frac{{{Df}}_{{\rm{R}}}}{d\lambda }=C[{f}_{{\rm{R}}}],\end{eqnarray} \tag{ 16 }$$

where λ is an affine parameter along a photon trajectory (geodesic) and $C[{f}_{{\rm{R}}}]$ accounts for interactions with matter: emission, absorption, and scattering of photons. The Liouville operator D/dλ is a derivative along the photon trajectory in phase space.

One can rewrite the Boltzmann equation as the radiative transfer equation,

$$\begin{eqnarray}&&\displaystyle \frac{D}{d\lambda }\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right)=\left(\displaystyle \frac{{\eta }_{\nu }}{{\nu }^{2}}\right)-\left(\nu {\chi }_{\nu }\right)\left(\displaystyle \frac{{I}_{\nu }}{{\nu }^{3}}\right).\end{eqnarray} \tag{ 17 }$$

Here the extinction coefficient

$$\begin{eqnarray}&&{\chi }_{\nu }\equiv {\alpha }_{\nu }+{\sigma }_{\nu },\end{eqnarray} \tag{ 18 }$$

and the emission coefficient

$$\begin{eqnarray}&&{\eta }_{\nu }\equiv {j}_{\nu }+{\eta }_{\nu }^{s}({I}_{\nu }),\end{eqnarray} \tag{ 19 }$$

where j_ν is the fluid emissivity, ${\eta }_{\nu }^{s}({I}_{\nu })$ is the scattering contribution to emissivity, σ_ν is the scattering extinction coefficient, and α_ν is the absorption coefficient. Each of the quantities in parentheses in (17) is invariant.

We neglect stimulated Compton scattering. The ratio of stimulated to spontaneous scattering is the photon occupation number in the scattered state. Models of highly sub-Eddington accretion onto supermassive black holes commonly feature: (1) relativistic electrons with Θ_e ≡ kT_e/(m_e c²) > 1, corresponding to a mean amplification factor after Compton scattering of $\approx 16{{\rm{\Theta }}}_{{\rm{e}}}^{2}$; and (2) a low frequency (millimeter or far-IR) peak in the spectrum at ν_pk where the synchrotron absorption optical depth is ${\mathcal{O}}(1)$. The energetically important single scattering events therefore produce scattered photons with ${\nu }_{\mathrm{sc}}\sim {\nu }_{\mathrm{pk}}16{{\rm{\Theta }}}_{{\rm{e}}}^{2}$. For moderate accretion rates (i.e., scattering depth τ_s < 1 for the disk), the photon occupation number at ν_sc is small, so stimulated Compton scattering will be negligible.

A consequence of our neglect of stimulated Compton scattering is that in a purely scattering medium the radiation field will approach a Wien (Boltzmann) distribution rather than a Bose–Einstein distribution. We verify this in Section 4.2.

To summarize, the governing equations for the radiation are (13), (14), and (17), together with appropriate expressions for the emission, scattering, and absorption coefficients.

In this section we adopt units such that c = 1 unless otherwise stated. It is apparent from Equations (3) and (17) that the fluid acts on the radiation through extinction and emission coefficients ${\chi }_{\nu }$ and η_ν. The radiation acts on the fluid through the four-force density G_μ. We want to make these representations consistent. Begin with the manifestly covariant expression

$$\begin{eqnarray}&&{R}^{\mu \nu }=\int \displaystyle \frac{{d}^{3}p}{\sqrt{-g}{p}^{t}}\;{p}^{\mu }{p}^{\nu }{f}_{{\rm{R}}}.\end{eqnarray} \tag{ 20 }$$

Then

$$\begin{eqnarray}&&{G}^{\mu }=-{{\rm{\nabla }}}_{\nu }\int \displaystyle \frac{{d}^{3}p}{\sqrt{-g}{p}^{t}}{p}^{\mu }{p}^{\nu }{f}_{{\rm{R}}}=-h\int \displaystyle \frac{{d}^{3}p}{\sqrt{-g}{p}^{t}}{p}^{\mu }\displaystyle \frac{{{Df}}_{{\rm{R}}}}{d\lambda }.\end{eqnarray} \tag{ 21 }$$

The last equality follows from an expansion of D/dλ and an integration by parts over momentum space (Lindquist 1966).

Using f_R = I_ν/(h⁴ν³), and Equations (17) and (21)

$$\begin{eqnarray}&&{G}^{\mu }=\displaystyle \frac{1}{{h}^{3}}\int \displaystyle \frac{{d}^{3}p}{\sqrt{-g}{p}^{t}}{p}^{\mu }\left[(\nu {\chi }_{\nu })({I}_{\nu }/{\nu }^{3})-({\eta }_{\nu }/{\nu }^{2})\right],\end{eqnarray} \tag{ 22 }$$

where ν, ${\chi }_{\nu }$, I_ν, and η_ν are all evaluated in a frame with four-velocity u^μ.

G^μ can be evaluated in an orthonormal tetrad frame comoving with the fluid

$$\begin{eqnarray}&&{e}_{(a)}^{\mu },\;\;{e}_{(t)}^{\mu }={u}^{\mu }.\end{eqnarray} \tag{ 23 }$$

We will call this the "fluid frame." In the fluid frame,

$$\begin{eqnarray}&&{G}_{(a)}=\int d\nu d{\rm{\Omega }}\left({\chi }_{\nu }{I}_{\nu }-{\eta }_{\nu }\right){n}_{(a)},\end{eqnarray} \tag{ 24 }$$

where ${n}_{(a)}\equiv {p}_{(a)}/(h\nu )$. Then

$$\begin{eqnarray}&&{G}^{\mu }={e}_{(a)}^{\mu }{G}^{(a)},\end{eqnarray} \tag{ 25 }$$

which is manifestly consistent with energy–momentum gains and losses by the radiation field.

The fluid integrator in bhlight is taken from harm, a conservative second order shock-capturing scheme on a two-dimensional mesh with an arbitrary spacetime metric. Here we give a brief summary of the method. Also, we adopt units such that c = 1, and for black holes we set GM = 1.

The fluid sector in bhlight updates a set of conserved variables ${\boldsymbol{U}}$:

$$\begin{eqnarray}&&{\boldsymbol{U}}=\sqrt{-g}\left({\rho }_{0}{u}^{t},{T}_{t}^{t},{T}_{i}^{t},{B}^{i}\right),\end{eqnarray} \tag{ 26 }$$

corresponding to the variables whose coordinate time derivatives are given in Section 2.1. These conserved variables are updated by fluxes ${\boldsymbol{F}}$:

$$\begin{eqnarray}&&{\boldsymbol{F}}=\sqrt{-g}\left({\rho }_{0}{u}^{i},{T}_{t}^{i},{T}_{j}^{i},{B}^{i}{\tilde{v}}^{j}-{B}^{j}{\tilde{v}}^{i}\right),\end{eqnarray} \tag{ 27 }$$

which in turn are calculated from the primitive variables ${\boldsymbol{P}}$:

$$\begin{eqnarray}&&{\boldsymbol{P}}=\left({\rho }_{0},u,{\tilde{v}}^{i},{B}^{i}\right),\end{eqnarray} \tag{ 28 }$$

where

$$\begin{eqnarray}&&{\tilde{v}}^{i}={v}^{i}+\displaystyle \frac{\gamma {\beta }^{i}}{\alpha },\end{eqnarray} \tag{ 29 }$$

where vⁱ = uⁱ/u⁰ is the fluid spatial three-velocity, $\gamma =\sqrt{1+{g}_{{ij}}{u}^{i}{u}^{j}}$, $\alpha =\sqrt{-1/{g}^{00}}$ is the lapse, and βⁱ = g⁰ⁱ α is the shift. Unlike vⁱ, ${\tilde{v}}^{i}$ ranges over $-\infty \lt {\tilde{v}}^{i}\lt \infty$. In the Newtonian formulation all transformations between the nonrelativistic analogs of ${\boldsymbol{U}}$, ${\boldsymbol{F}}$, and ${\boldsymbol{P}}$ are analytic, but in the covariant formulation there is no general analytic form for ${\boldsymbol{P}}$ $({\boldsymbol{U}})$.

The fluid update each timestep maps ${\boldsymbol{P}}$ⁿ to its updated value ${\boldsymbol{P}}$^{n + 1} by updating the conserved variables. Beginning with ${\boldsymbol{P}}$ⁿ, the scheme calculates ${\boldsymbol{U}}$ⁿ = ${\boldsymbol{U}}$ (${\boldsymbol{P}}$ⁿ) and ${\boldsymbol{F}}$ⁿ = ${\boldsymbol{F}}$ (${\boldsymbol{P}}$ⁿ) via closed-form expressions, for ${\boldsymbol{F}}$ⁿ after a reconstruction step that estimates ${\boldsymbol{P}}$ⁿ at zone boundaries from values at zone centers. The update ${{\boldsymbol{U}}}^{n}\to {{\boldsymbol{U}}}^{n+1}$ over a timestep Δt is given by

$$\begin{eqnarray}{{\boldsymbol{U}}}^{n+1} & = & {{\boldsymbol{U}}}^{n}+{\rm{\Delta }}t\left(-\displaystyle \frac{{{\boldsymbol{F}}}_{i+1,j}^{n+1/2}-{{\boldsymbol{F}}}_{i,j}^{n+1/2}}{{\rm{\Delta }}{x}^{1}}\right.\\ & & -\left.\displaystyle \frac{{{\boldsymbol{F}}}_{i,j+1}^{n+1/2}-{{\boldsymbol{F}}}_{i,j}^{n+1/2}}{{\rm{\Delta }}{x}^{2}}+{\dot{{\boldsymbol{U}}}}^{n+1/2}\right),\end{eqnarray} \tag{ 30 }$$

where $\dot{{\boldsymbol{U}}}$ represents the source terms such as those associated with the spacetime connection, values at n + 1/2 are estimated from a similar first-order step to ${\boldsymbol{U}}$^{n + 1/2}, and i, j here denote spatial indices in x¹ and x², respectively. This forms a second order, explicit timestepping scheme to t + Δ t, and then ${\boldsymbol{P}}$^{n + 1} is found by numerically solving ${\boldsymbol{U}}$ (${\boldsymbol{P}}$^{n + 1}) = ${\boldsymbol{U}}$^{n + 1} (see Noble et al. 2006 and Mignone & McKinney 2007).

The fluxes ${\boldsymbol{F}}$(${\boldsymbol{P}}$) are evaluated at zone faces using Local Lax Friedrichs fluxes. Primitive variables on either side of the zone face are determined through slope-limited linear reconstruction. We typically use the monotonized central limiter for reconstruction, but it is trivial to use higher order methods as well.

Naively differencing the induction equation (10) will not preserve a numerical representation of the no-monopoles constraint (11); the monopole density will undergo a random walk from zero with a step size determined by truncation error. Unless directly controlled, the monopole density can grow quickly and corrupt the solution. A variety of techniques for avoiding or removing magnetic monopoles exist; bhlight employs the flux-interpolated constrained transport (flux-CT) scheme introduced by Tóth (2000). Although this introduces some additional diffusivity into the scheme, it is simple and effective. Details of the implementation are given in Gammie et al. (2003).

bhlight uses nearly the same Monte Carlo implementation as grmonty, with a few important differences. The Monte Carlo samples are referred to here as superphotons. Each superphoton has a weight w (the number of photons carried by the superphoton), a momentum p_μ (${p}_{\mu }={\rm{\hslash }}{k}_{\mu }$ and p_μp^μ = 0), and a position x^μ.

The Monte Carlo representation of the the photon distribution function is

$$\begin{eqnarray}&&{f}_{{\rm{R}},\mathrm{MC}}=\underset{k}{\displaystyle \sum }\;{w}_{k}\;{\delta }^{3}({x}^{i}-{x}_{k}^{i}){\delta }^{3}({p}_{j}-{p}_{j,k}),\end{eqnarray} \tag{ 31 }$$

where ${\delta }^{3}({x}^{i}-{x}_{k}^{i})=\delta ({x}^{1}-{x}_{k}^{1})\delta ({x}^{2}-{x}_{k}^{2})\delta ({x}^{3}-{x}_{k}^{3})$, etc. The sum is taken over all photon samples in the model, labeled by the index k, and w_k are the weights. Like f_R, f_R,MC is invariant because w_k, ${\delta }^{3}({x}^{i}-{x}_{k}^{i})/(\sqrt{-g}{p}^{t})$, and ${\delta }^{3}({p}_{j}-{p}_{j,k})\sqrt{-g}{p}^{t}$ (with p_j covariant) are all invariant.

Using Equation (20), the stress–energy tensor is

$$\begin{eqnarray}&&{R}^{\mu \nu }=\underset{k}{\displaystyle \sum }\displaystyle \frac{{p}_{k}^{\mu }{p}_{k}^{\nu }}{\sqrt{-g}{p}_{k}^{t}}{w}_{k}{\delta }^{3}({x}^{i}-{x}_{k}^{i}).\end{eqnarray} \tag{ 32 }$$

This can be averaged over a three-volume Δ³x = Δx¹ Δx² Δx³ to obtain an estimate for ${\bar{R}}^{\mu \nu }$:

$$\begin{eqnarray}&&{\bar{R}}^{\mu \nu }\approx \displaystyle \frac{1}{{{\rm{\Delta }}}^{3}x}\int {d}^{3}{{xR}}^{\mu \nu }=\displaystyle \frac{1}{\sqrt{-g}{{\rm{\Delta }}}^{3}x}\underset{k}{\displaystyle \sum }\displaystyle \frac{{p}_{k}^{\mu }{p}_{k}^{\nu }}{{p}_{k}^{t}}{w}_{k}\end{eqnarray} \tag{ 33 }$$

where now the sum is taken only over photons within the three-volume (zone) in question.

3.2.1. Initializing the Radiation Field

How should one initialize f_R,MC? In bhlight's target applications this question usually does not arise because f_R relaxes rapidly to a quasi-equilibrium, so one can set f_R = 0 in the initial conditions. In test problems, however, an accurate initial f_R may be required. In this case one wants to sample a set of photons in a single zone centered at x_c; that is, we want to sample Δ³xf_R(x_c).

One strategy is to sample f_R directly in a coordinate frame, using the invariance of f_R. For example, if f_R is thermal in the fluid frame, then the distribution function in any coordinate frame is ${c}^{2}{B}_{\nu }/({h}^{4}{\nu }^{3})$, where B_ν is the Planck function, and $\nu =-{u}^{\mu }{k}_{\mu }/(2\pi )$.

A second strategy, which we adopt, is to sample f_R in the fluid frame (comoving indices are denoted by parentheses). Then we must take care: ${{\rm{\Delta }}}^{3}{{xf}}_{{\rm{R}}}({x}_{c})$ is not invariant, because the volume element Δ³x is not invariant. But ${{\rm{\Delta }}}^{3}x\sqrt{-g}{p}^{t}$ is invariant, so $({{\rm{\Delta }}}^{3}x/{{\rm{\Delta }}}^{3}x\prime )=\sqrt{-g}{p}^{t}/{p}^{(t)},$ where Δ³x' is a fluid frame volume element, and $\sqrt{-g}=1$ in the fluid frame. Then ${{\rm{\Delta }}}^{3}{{xf}}_{{\rm{R}}}({x}_{c})=({{\rm{\Delta }}}^{3}x/{{\rm{\Delta }}}^{3}x\prime ){{\rm{\Delta }}}^{3}x\prime {f}_{{\rm{R}}}=(\sqrt{-g}{p}^{t}/{p}^{(t)}){{\rm{\Delta }}}^{3}x\prime {f}_{{\rm{R}}}$. This suggests that we can sample Δ³x' f_R in the fluid frame and then multiply the photon number dN by the corresponding $\sqrt{-g}{p}^{t}/{p}^{(t)}$ to obtain a fair sample of Δ³xf_R in the coordinate frame.

This second strategy can be described more explicitly in terms of the Monte Carlo samples as follows. A list of photons in a single zone is obtained by taking

$$\begin{eqnarray}&&{\displaystyle \int }_{{{\rm{\Delta }}}^{3}x}\;{d}^{3}x\;{f}_{{\rm{R}},\mathrm{MC}}=\underset{k}{\displaystyle \sum }\;{w}_{k}\;{\delta }^{3}({p}_{j}-{p}_{j,k}),\end{eqnarray} \tag{ 34 }$$

where the sum is over photon samples in a single zone. This is not invariant because δ³(p_j − p_j,k) is not invariant, so one must take care in sampling p_j,k and w_k. Suppose we sample f_RΔ³x' in the tetrad frame. This gives us a list of weights and momenta. We can transform back to the fluid frame using the invariance of $\sqrt{-g}{p}^{t}\delta ({p}_{j}-{p}_{j,k}),$ so that each w_kδ(p_j − p_j,k) in the tetrad frame becomes ${w}_{k}\sqrt{-g}({p}^{t}/{p}^{(t)})\delta ({p}_{j}-{p}_{j,k})$ in the coordinate frame. We can therefore obtain a fair sample by adjusting the weights by a factor of $\sqrt{-g}{p}^{t}/{p}^{(t)}$ in transforming from the fluid frame to the coordinate frame.

3.2.2. Geodesic Integration

3.2.3. Units

The radiation sector uses cgs units, except that photon wavevector components are measured in units of the electron rest mass energy. We therefore have to convert between units in the fluid sector and units in the radiation sector. For black hole spacetimes, this is done by choosing a cgs value for the fluid length unit and the fluid time unit, here

$$\begin{eqnarray}&&{\mathcal{L}}=\displaystyle \frac{{GM}}{{c}^{2}},\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&{\mathcal{T}}=\displaystyle \frac{{GM}}{{c}^{3}},\end{eqnarray} \tag{ 36 }$$

respectively. We also need a mass unit. Notice that the mass unit is not provided by the black hole mass in the test fluid $(\rho {{\mathcal{L}}}^{3}\ll M)$ approximation used here. Instead we must scale the density, or equivalently the mass accretion rate, by choosing a cgs value for the mass unit ${\mathcal{M}}.$ Then, e.g., ${\rho }_{\mathrm{CGS}}={\rho }_{\mathrm{FLUID}}{\mathcal{M}}/{{\mathcal{L}}}^{3}.$

The components of the code photon wavevector k^μ are measured in units of m_ec². One might then be concerned about consistency between the transfer equation and the geodesic equation. The only condition for consistency is that the differential optical depth dτ_ν = (νκ_ν) dλ, which in turn requires that the νdλ = ds, i.e., that the units used in defining ν and dλ be consistent and that the correct conversion be made from fluid sector units to cgs. In practice, then, we evaluate νκ_ν in cgs units in the fluid frame and set $d{\tau }_{\nu }={\mathcal{N}}(\nu {\kappa }_{\nu })d\lambda$, with ${\mathcal{N}}\equiv h{\mathcal{L}}/({m}_{{\rm{e}}}{c}^{2})$.

3.2.4. Superphoton Weighting

3.2.5. Emissivity

At each timestep we sample the emissivity in the invariant four-volume $\sqrt{-g}{\rm{\Delta }}t{{\rm{\Delta }}}^{3}x$ of each zone based on the fluid values at the half-step. It is easiest to sample the fluid emission in a comoving tetrad, where we have a simple expression for the emissivity.

The emissivity is

$$\begin{eqnarray}&&{j}_{\nu }=\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{{dE}}{{d}^{3}{xdtd}\nu d{\rm{\Omega }}}.\end{eqnarray} \tag{ 37 }$$

Using ${dE}=h\nu {dN}$ where N is the number of photons, we can write

$$\begin{eqnarray}&&\int d{\rm{\Omega }}\;{j}_{\nu }=\displaystyle \frac{h}{\sqrt{-g}}\displaystyle \frac{{dN}}{{d}^{3}{xdtd}\mathrm{log}\nu }.\end{eqnarray} \tag{ 38 }$$

Since dN = w dN_s where N_s is the number of superphotons, we can then write for the number of superphotons produced per logarithmic interval in a single zone with volume Δ³x:

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{{\rm{s}}}}{d\mathrm{log}\nu }=\sqrt{-g}{\rm{\Delta }}t{{\rm{\Delta }}}^{3}x\displaystyle \frac{1}{{hw}(\nu )}\int d{\rm{\Omega }}\;{j}_{\nu }.\end{eqnarray} \tag{ 39 }$$

In writing (39) we have made use of the invariance of $\sqrt{-g}{d}^{3}{xdt}$. Here w(ν) ∝ 1/ν, and the constant of proportionality is set dynamically to keep the number of superphotons in the computational domain approximately constant.

bhlight samples Equation (39) between minimum and maximum frequencies ν_min and ν_max, where the limits are set so that νj_ν is small outside this region in frequency space. It then uses a rejection scheme, sampling a uniform distribution in $\mathrm{log}{\nu }_{\mathrm{min}}\lt \mathrm{log}\nu \lt \mathrm{log}{\nu }_{\mathrm{max}}$ and a uniform distribution in $0\lt r\lt {({{dN}}_{{\rm{s}}}/d\mathrm{log}\nu )}_{\mathrm{max}}$, where ${({{dN}}_{{\rm{s}}}/d\mathrm{log}\nu )}_{\mathrm{max}}$ is the maximum of Equation (39). A sample is rejected when $r\gt {{dN}}_{{\rm{s}}}/d\mathrm{log}\nu$.

The angular distribution of photons is also sampled by rejection. The photon direction is given by (θ,ϕ), where θ is the angle between the magnetic field and photon direction in the fluid frame and ϕ is the corresponding azimuthal angle. bhlight samples a uniform distribution in $0\leqslant \mathrm{cos}\theta \lt 1$, and a uniform distribution in 0 ≤ r < 1. A sample is rejected if r > j_ν(θ)/j_ν,max. It then samples a uniform distribution in 0 ≤ ϕ < 2π. To ensure that the net force due to emission in the fluid frame is zero to machine precision, photons are generated in pairs. Thus, a second photon is generated with the same frequency and $\mathrm{cos}{\theta }^{\prime }=-\mathrm{cos}\theta$ and ϕ' = ϕ + π. In the fluid frame, k^t = ω, ${k}^{(x)}=\omega \mathrm{sin}\theta \mathrm{cos}\phi$, ${k}^{(y)}=\omega \mathrm{sin}\theta \mathrm{sin}\phi$, and ${k}^{(z)}=\omega \mathrm{cos}\theta$, where ${e}^{(z)}$ is parallel to the magnetic field. Once we have a superphoton sample in the comoving frame it is transformed to the coordinate frame using a pre-constructed orthonormal tetrad ${e}_{(a)}^{\mu }$. The superphoton x^μ is set to the zone center to avoid additional orthonormal tetrad construction.⁶

Sampling is a subdominant computational expense in bhlight, so although one could develop more efficient sampling schemes, a simple rejection scheme is adequate.

Four-momentum is locally conserved and so superphoton emission implies a back-reaction on the emitting fluid. A pair of emitted superphotons with wavevectors ${k}_{1}^{\mu }$, ${k}_{2}^{\mu }$ correspond to a change in four-momentum Δp^μ (in fluid code units):

$$\begin{eqnarray}&&{\rm{\Delta }}{p}^{\mu }={\mathcal{P}}w\left({k}_{1}^{\mu }+{k}_{2}^{\mu }\right),\end{eqnarray} \tag{ 40 }$$

where ${\mathcal{P}}={m}_{{\rm{e}}}/{\mathcal{M}}$, which in turn specifies the contribution to the four-force density ΔG_μ:

$$\begin{eqnarray}&&{\rm{\Delta }}{G}_{\mu }=-\displaystyle \frac{1}{\sqrt{-g}{{\rm{\Delta }}}^{3}x{\rm{\Delta }}t}{\rm{\Delta }}{p}_{\mu },\end{eqnarray} \tag{ 41 }$$

where all geometric quantities are evaluated at zone centers. Because photons are emitted in pairs, the spatial components of k^μ cancel in the fluid frame, and δp_μ ∝ u_μ.

3.2.6. Absorption

grmonty treats absorption deterministically by continuously decrementing w along a ray. A similar deterministic procedure has been shown to suppress noise in Monte Carlo radiation hydro schemes (e.g., Noebauer et al. 2012) in flat space. However, formulating such a scheme in general relativity, where photons move along geodesics, is more complicated because the photons follow curved trajectories through each zone.

In bhlight we treat absorption probabilistically. While stepping a superphoton by Δλ along a geodesic, the incremental optical depth to absorption ${\rm{\Delta }}{\tau }_{a}={\mathcal{N}}{\kappa }_{\nu ,\mathrm{abs}}\nu {\rm{\Delta }}\lambda$. Here ν κ_ν,abs is the invariant absorption coefficient, evaluated in the fluid frame. An absorption occurs if

$$\begin{eqnarray}&&{\rm{\Delta }}{\tau }_{a}\gt -\mathrm{log}{r}_{a},\end{eqnarray} \tag{ 42 }$$

where 0 < r_a < 1 is sampled uniformly; the absorption occurs at ${\rm{\Delta }}{\tau }_{a}=-\mathrm{log}{r}_{a}$. To process the event we push the superphoton back, $\lambda \to \lambda +{\rm{\Delta }}\lambda (\mathrm{log}{r}_{a}/{\rm{\Delta }}{\tau }_{a})$, put the superphoton four-momentum into the fluid at that location, and annihilate the superphoton.

The four-momentum change in the fluid Δp^μ due to absorption of a superphoton with wavevector k^μ is

$$\begin{eqnarray}&&{\rm{\Delta }}{p}^{\mu }={\mathcal{P}}{{wk}}^{\mu }.\end{eqnarray} \tag{ 43 }$$

This can be expressed as a contribution to the radiation four-force density ΔG_μ:

$$\begin{eqnarray}&&{\rm{\Delta }}{G}_{\mu }=\displaystyle \frac{{\rm{\Delta }}{p}_{\mu }}{\sqrt{-g}{{\rm{\Delta }}}^{3}x{\rm{\Delta }}t},\end{eqnarray} \tag{ 44 }$$

where $\sqrt{-g}$ is evaluated at the zone center.

3.2.7. Scattering

We treat scattering probabilistically in bhlight, as in grmonty. Scattering is similar to absorption, i.e., scattering occurs when

$$\begin{eqnarray}&&{\rm{\Delta }}{\tau }_{{\rm{s}}}\gt -(\mathrm{log}{r}_{{\rm{s}}})/{b}_{{\rm{s}}},\end{eqnarray} \tag{ 45 }$$

where 0 < r_s < 1 is sampled uniformly, ${\rm{\Delta }}{\tau }_{{\rm{s}}}={\mathcal{N}}{\kappa }_{\nu ,{\rm{s}}}\nu {\rm{\Delta }}\lambda$, and νκ_ν,s is the invariant scattering opacity. Because scattering events are rare but energetically important we have introduced a bias parameter b_s > 1 to enhance the probability of sampling scattering events. To process the event, we push the photon back along the geodesic from λ + Δλ to $\lambda +{\rm{\Delta }}\lambda \mathrm{log}{r}_{{\rm{s}}}/({b}_{{\rm{s}}}{\rm{\Delta }}{\tau }_{{\rm{s}}})$. To preserve photon number a scattered superphoton is created with weight w_s = w/b_s and the original superphoton has weight set to w' = w − w_s.

In general, a superphoton is subject to both absorption and scattering simultaneously. In a deterministic treatment, the code must dynamically choose which process, if any, to apply to the superphoton. To handle this in an unbiased manner, for each photon assuming that at least one of the inequalities Equations (42) and (45) has been satisfied, we choose which interaction to process according to a similar weighted sampling. That is, for

$$\begin{eqnarray}&&\displaystyle \frac{-\mathrm{log}{r}_{a}}{{\rm{\Delta }}{\tau }_{a}}\lt \displaystyle \frac{-\mathrm{log}{r}_{{\rm{s}}}}{{b}_{{\rm{s}}}{\rm{\Delta }}{\tau }_{{\rm{s}}}},\end{eqnarray} \tag{ 46 }$$

the absorption interaction is chosen; elsewhere, scattering is chosen. With this method, a large optical depth or bias in one interaction will not serve to decrease the physical effect of the other interaction, although it will increase the number of superphotons required to resolve both interactions simultaneously.

How should we set b_s? Most of our models have ${\tau }_{{\rm{s}}}\ll 1$, so only ≈τ_s superphotons would produce scattering events if b_s = 1, and the energy per superphoton would increase by the mean amplification factor⁷ $A\approx 1+4{{\rm{\Theta }}}_{{\rm{e}}}+16{{\rm{\Theta }}}_{{\rm{e}}}^{2}$. This suggests that we should set b_s ∼ A to maintain constant energy per superphoton. There are two failure modes to be avoided, however. First, if b_s ≳ 1/τ_s then each superphoton will scatter more than once and the number of superphotons on the grid will grow exponentially. Second, if Awhν is larger than the total internal energy in the zone where the scattering occurs then the zone energy will be negative after scattering (more on this in the next subsection). Together, this suggests that we set

$$\begin{eqnarray}&&{b}_{{\rm{s}}}=A\mathrm{MAX}\left({\mathcal{C}},\displaystyle \frac{{wh}\nu }{u\sqrt{-g}{{\rm{\Delta }}}^{3}x}\right).\end{eqnarray} \tag{ 47 }$$

Here ${\mathcal{C}}$ depends only on t and is dynamically adjusted to control the number of scattered superphotons in the simulation. The requirement b_sτ_s ∼ 1 is equivalent to each superphoton scattering once between emission and escape through the boundary (neglecting absorption). Over a timestep Δt, one can estimate the number of photons that escape through the boundaries of a domain with linear dimension L as NcΔt/L, where N is the desired total number of superphotons (which sets the weights for emission as described previously). To enforce the requirement that each superphoton scatter once, ${\mathcal{C}}$ is calculated dynamically as the ratio of this estimate to the real number of scattering events per timestep, averaged over some timescale.

Each scattered superphoton is generated from an incident superphoton wavevector k^μ as follows. The four-momentum p^μ of the scattering electron is sampled from a thermal (Maxwell–Jüttner) distribution according to the procedure described in Canfield et al. (1987). The scattered superphoton wavevector ${k}_{{\rm{s}}}^{\mu }$ is sampled from the Klein–Nishina differential scattering cross section in the rest frame of the scattering electron and boosted to the fluid frame and then transformed to the coordinate frame. It is then assigned a weight and entered in the list of active superphotons.

Each scattering event generates a change in fluid four-momentum,

$$\begin{eqnarray}&&{\rm{\Delta }}{p}^{\mu }={\mathcal{P}}\displaystyle \frac{w}{{b}_{{\rm{s}}}}\left({k}^{\mu }-{k}_{{\rm{s}}}^{\mu }\right)\end{eqnarray} \tag{ 48 }$$

and a corresponding contribution to the fluid through the four-force density G_μ,

$$\begin{eqnarray}&&{\rm{\Delta }}{G}_{\mu }=\displaystyle \frac{{\rm{\Delta }}{p}_{\mu }}{\sqrt{-g}{{\rm{\Delta }}}^{3}x{\rm{\Delta }}t},\end{eqnarray} \tag{ 49 }$$

where $\sqrt{-g}$ is evaluated at the zone center.

The radiation force is treated in an operator-split fashion. The fluid integrator initially updates the conserved variables ${\boldsymbol{U}}$ from step n to $n+1$ over the entire timestep Δt without radiation, i.e., it performs ${{\boldsymbol{U}}}^{n}\to {{\boldsymbol{U}}}^{n+1\prime }$ as described in Section 3.1. This fluid integration generates half-step fluid primitive variables ${{\boldsymbol{P}}}^{n+1/2}$; these values are sent to the radiation sector and used to evaluate the total radiation four-force density ${G}_{\mu }({{\boldsymbol{P}}}^{n+1/2})$ for each zone. The fluid integrator then updates the fluid variables with the radiation interaction ${{\boldsymbol{U}}}^{n+1\prime }\to {{\boldsymbol{U}}}^{n+1}$ by considering only the radiation contribution to the conserved energy and momentum variables:

$$\begin{eqnarray}&&{\left(\sqrt{-g}{T}_{\nu }^{t}\right)}^{n+1}={\left(\sqrt{-g}{T}_{\nu }^{t}\right)}^{n+1\prime }+{\rm{\Delta }}t\sqrt{-g}{G}_{\nu }.\end{eqnarray} \tag{ 50 }$$

${\boldsymbol{U}}$ⁿ⁺¹ then are the final conserved fluid variables at the $(n+1)\mathrm{th}$ step. The evolution is therefore first order in time.

The evolution is explicit and the radiation and fluid share a common timestep Δt, which we set to the minimum grid zone light crossing time ≈Δx/c, where Δx is a characteristic zone lengthscale. As is well known, the radiation source terms are stiff when the timescale for exchange of energy–momentum between the fluid and radiation is smaller than a timestep. The cooling time τ_cool ≡ u/Λ, where Λ ≡ cooling rate per unit volume, $=\;{u}_{\mu }{G}^{\mu }\sim {u}_{\mathrm{rad}}c/{\lambda }_{\mathrm{mfp}}$, where λ_mfp is a suitably frequency-averaged absorption mean free path and ${u}_{\mathrm{rad}}={R}^{\mu \nu }{u}_{\mu }{u}_{\nu }$ is the radiation energy density in the fluid frame (one can perform a similar estimate for Compton cooling). Thus the source term is stiff if u/Λ < Δx/c or $({u}_{\mathrm{rad}}/u)({\rm{\Delta }}x/{\lambda }_{\mathrm{mfp}})\gt 1$, or when the optical depth across a zone exceeds u/u_rad.

For our scheme we must also consider robustness in the presence of Monte Carlo noise. Even if τ_cool/Δt > 1 the cooling rate may fluctuate upward so that a zone loses all its thermal energy in a single timestep. This can happen if $u/{\langle {{\rm{\Lambda }}}^{2}\rangle }^{1/2}\lesssim {\rm{\Delta }}x/c$. This will differ from the usual stiffness condition only when the number of absorption events per timestep in a zone is small compared to one. The condition for robustness against this failure mode, "supercooling," where a single photon causes the zone to lose all its internal energy, is that ${wh}\nu \lt u{{\rm{\Delta }}}^{3}x$ (where we have left out geometric factors).

Where, then, will bhlight fail? The radiation force source terms are stiff when the optical depth across a single zone exceeds u/u_rad. For black hole accretion applications we expect this only for models in the high accretion rate regime, $\dot{M}\gtrsim {10}^{-3}{\dot{M}}_{\mathrm{Edd}}$, although the precise condition will depend on details of the evolution and the numerical setup. This problem could be remedied by using an implicit update, but Monte Carlo is probably not the optimal method for studying this regime anyway. The supercooling problem is more relevant for our target application to intermediate accretion rate black holes, and arises if the internal energy content of a zone is small compared to the typical superphoton energy. This can occur in low density regions over the poles of the black hole, but the fluid evolution is inaccurate there in any case (because the truncation error in internal energy is dominated by the magnetic field evolution, to which it is coupled via the total energy density), and negative internal energies are dealt with by harm's floor routines, resulting in a small nonconservation of energy.

bhlight is a hybrid MPI/OpenMP scheme in which a single node handles the fluid integration, multiple nodes evolve the radiation, and a single additional node acts as gatekeeper between the fluid and radiation sectors. During each timestep, the only exchanges are an array of radiation four-force densities to the fluid sector and an array of fluid variables to the radiation sector via the gatekeeper node. The gatekeeper node distributes the fluid variables to all radiation nodes, and reduces the four-force density contributions from each radiation node.

We evolve radiation on each node independently from other radiation nodes. After globally scaling the emission weights and scattering bias to yield approximately the desired number of superphotons at saturation, the code samples emission events on each radiation node, each of which has access to the entire array of fluid variables and maintains its own set of superphoton samples. Emission, absorption, and scattering events generate a four-force density contribution. At the end of every timestep, these contributions are reduced by the primary radiation node over MPI. Each radiation node is individually parallelized under OpenMP, further dividing the superphoton calculations across individual compute cores. We parallelize the main compute loops for the fluid sector with OpenMP, which enables completion in a reasonable clock time for an axisymmetric calculation.

In attempting to describe our numerical implementation in a coherent narrative we have omitted certain secondary topics, which we now collect here.

• 1.  
  The radiation sector in bhlight makes extensive use of random numbers. We use the Mersenne Twister algorithm from the GNU Scientific Library, with a different random seed for each MPI node and each OpenMP thread.
• 2.  
  In axisymmetric disk calculations, we implement a form of static mesh refinement by using modified Kerr–Schild (MKS) coordinates {t, x¹, x²}. x¹ and x² are related to the Kerr–Schild r, θ by $r=\mathrm{exp}\left({x}^{1}\right)+{r}_{0}$ and $\theta =\pi {x}^{2}+((1-{h}_{{\rm{s}}})/2)\mathrm{sin}(2\pi {x}^{2})$, where ${r}_{0}\in [0,\infty )$ and ${h}_{{\rm{s}}}\in (0,1]$ are free parameters.
• 3.  
  Truncation error in the geodesic integration causes k^μ to drift off the lightcone (this is a consequence of our decision to integrate all four components). We destroy superphotons with negative frequency in the fluid frame; for torus runs as in Section 5.2, we find ∼1.1 × 10⁻⁶ destructions per geodesic update. This problem does not occur in Cartesian coordinates in Minkowski space.
• 4.  
  Scattered superphotons may scatter any number of additional times during the same timestep, provided sufficient optical depth to do so.
• 5.  
  bhlight does not conserve momentum and energy to machine precision because of truncation error in the geodesic integrator. However, for the integrator tolerance and typical superphoton resolutions this is not significant (for torus runs as in Section 5.2). At some time, the average fractional error in energy relative to the initial energy at emission is ∼few × 10⁻⁷; if this were to become a leading source of error, increasing the integrator tolerance is not a significant expense.

We consider the temperature evolution of an optically thin, radiating, stationary, ideal, and homogeneous gas initially at temperature T₀. The gas obeys a γ-law equation of state i.e., p = (γ − 1)u. The density and velocity of the gas are fixed; only the temperature is allowed to evolve. The bremsstrahlung-like emissivity is

$$\begin{eqnarray}&&{j}_{\nu }={{Nn}}^{2}{T}^{-1/2}\mathrm{exp}\left(-h\nu /{k}_{{\rm{B}}}T\right),\end{eqnarray} \tag{ 51 }$$

where $N=5.4\times {10}^{-39}\;{\mathrm{cm}}^{3}\;{{\rm{K}}}^{1/2}\;{{\rm{s}}}^{-1}\;{\mathrm{Sr}}^{-1}\;{\mathrm{Hz}}^{-1}$ is a constant. The associated cooling rate Λ is given by

$$\begin{eqnarray}&&{\rm{\Lambda }}=-\displaystyle \frac{{du}}{{dt}}=\int d\nu d{\rm{\Omega }}{j}_{\nu }=4\pi \displaystyle \frac{{k}_{{\rm{B}}}{n}^{2}N}{h}{T}^{1/2},\end{eqnarray} \tag{ 52 }$$

which implies a temperature evolution

$$\begin{eqnarray}&&T(t)={T}_{0}{\left(1-\displaystyle \frac{t}{{t}_{f}}\right)}^{2}\end{eqnarray} \tag{ 53 }$$

valid from t_i = 0 to ${t}_{f}={{hT}}_{0}^{1/2}/((\gamma -1)\pi {Nn})$, the time at which the temperature of the fluid reaches zero.

For this realization we choose γ = 5/3, t_f = 10⁸ s, and T₀ = 10⁸ K. Figure 1 shows the resulting numerical evolution plotted against the analytic solution. Convergence in the L¹ norm,

$$\begin{eqnarray}&&{L}^{1}(f)\equiv \underset{i}{\displaystyle \sum }\left|{f}_{\mathrm{num},i}-{f}_{\mathrm{reference},i}\right|\end{eqnarray} \tag{ 54 }$$

is expected to scale as ${N}_{{\rm{s}}}^{-1/2}$; this can be seen in Figure 2, which is evaluated near t = t_f.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consider a closed, one-zone model in which Compton scattering is the only permitted interaction between an ideal, γ = 5/3 gas initially at temperature T_g,i and a swarm of photons all with initial frequency ν = ν₀. Fluid motion is suppressed; only the internal energy is allowed to evolve. The number of photons is conserved, and in thermal equilibrium T_g,f = T_r,f = T_f, the radiation approaches the Wien distribution,

$$\begin{eqnarray}&&{f}_{{\rm{R}}}=\displaystyle \frac{{n}_{\gamma }}{8\pi }{\left(\displaystyle \frac{c}{{{kT}}_{f}}\right)}^{3}\mathrm{exp}(-h\nu /({k}_{{\rm{B}}}{T}_{f})),\end{eqnarray} \tag{ 55 }$$

where n_γ is the number density of photons. This tests the scattering kernel and Compton heating and cooling of the fluid.

We set the electron number density n = 2.5 × 10¹⁷ cm⁻³, T_g,i = 5 × 10⁷ K, ν₀ = 3 × 10¹⁶ Hz, and n_γ = 2.38 × 10¹⁸ cm⁻³. The characteristic (Compton) relaxation time is the photon mean free time ${({n}_{{\rm{e}}}{\sigma }_{T}c)}^{-1}$ divided by the fractional energy change per scattering, $\sim {kT}/({m}_{{\rm{e}}}{c}^{2})$, yielding a Comptonization timescale ≃0.02 s. We can predict the final state using: (1) thermal equilibrium; (2) conservation of photon number; (3) that the final photon distribution is Boltzmann; (4) conservation of total energy. This yields T_f = 5.19 × 10⁶ K. Figure 3 shows fluid and gas temperatures equilibrating at the correct temperature on approximately the estimated timescale, and ${u}_{\nu }\equiv {dE}/({d}^{3}x\;d\mathrm{log}\nu )$ plotted against the anticipated distribution along with associated residuals. Note that we define ${T}_{r}\equiv {E}_{r}/(3{n}_{\gamma }{k}_{{\rm{B}}})$ (E_r is the radiation energy density), which assumes a Wien distribution, and so this value is strictly valid only at late time.

Zoom In Zoom Out Reset image size

As a test of the full transfer equation in one-dimensional (1D) with gray absorption, consider a sinusoidal temperature perturbation in a static gas. Mihalas & Mihalas (1984) show that in this case the full transfer equation plus gas energy equation admit damped solutions. The eigenmode has perturbations $\propto \mathrm{exp}(i({kx}-t/{t}_{\mathrm{RR}}))$ for wavenumber k and decay time t_RR. The dispersion relation is

$$\begin{eqnarray}&&{t}_{\mathrm{RR}}(k)={\left[\displaystyle \frac{16\sigma \kappa {T}_{0}^{3}}{\rho {c}_{v}}\left(1-\displaystyle \frac{{\alpha }_{0}}{k}{\mathrm{cot}}^{-1}\left(\displaystyle \frac{{\alpha }_{0}}{k}\right)\right)\right]}^{-1},\end{eqnarray} \tag{ 56 }$$

where T₀ is the mean temperature, ρ is the material density, α₀ is the frequency-integrated extinction coefficient, and c_v is the specific heat capacity. We simulate this problem in bhlight by initializing one wavelength in a 1D box in local radiative equilibrium with 64 grid zones and periodic boundary conditions. The amplitude of the initial perturbation is 0.05 T₀. We evolve this system for a variety of optical depths per wavelength τ by varying α₀. We obtain a decay time from the amplitude of the best fit sinusoid at t = t_RR. We find good agreement with Equation (56) in both optically thin and optically thick regimes, as shown in Figure 4.

Zoom In Zoom Out Reset image size

We also examine convergence for τ = 1 (with 100 grid zones); this result is shown in Figure 5. Evidently the errors scale as ${N}_{{\rm{s}}}^{-1/2}$, as expected.

Zoom In Zoom Out Reset image size

We now also consider linear modes of the full equations of 1D radiation magnetohydrodynamics; that is, we now include momentum exchange and magnetic fields. Acquiring even a linear solution to these equations with full transport is challenging, and so we resort to the approximate, relativistic Eddington closure scheme of Farris et al. (2008) with gray opacity κ from which to extract plane-wave solutions ${\boldsymbol{P}}$. Details of our calculation for a perturbation δ${\boldsymbol{P}}$ for $\delta \sim \mathrm{exp}(\omega t-{ikx})$ about a thermal equilibrium ${\boldsymbol{P}}$₀ are given in Appendix C. With a magnetic field ${B}^{i}=({B}_{0},{B}_{0},0)$ we follow the nonrelativistic treatment of Jiang et al. (2012) by confining variation to a plane (suppressing Alfvèn waves): ${\boldsymbol{P}}=(\rho ,u,{u}^{1},{u}^{2},{B}^{1},{B}^{2},E,{F}^{1},{F}^{2})=({\rho }_{0}+\;\delta \rho ,{u}_{0}+\delta u,$ $\delta {u}^{1},\delta {u}^{2},{B}_{0},{B}_{0}+\delta {B}^{2},{E}_{0}+\delta E,\delta {F}^{1},\delta {F}^{2})$. bhlight is not designed to evolve perturbed equilibria; we focus on cases which accomodate both bhlight's numerical limitations as well as the discrepancy between the Eddington closure and bhlight's full transport (i.e., we consider only many optical depths per wavelength). We study convergence of two specific cases with significant radiation pressure: a nonrelativistic radiation-modified slow MHD mode, and a relativistic radiation-modified fast MHD mode. In this section, c = k_B = 1. For all calculations, we use a box of length L = 1 with 128 grid zones, evolve to final time t_f, set the wavenumber of the perturbation k = 2 π, and normalize the δ${\boldsymbol{P}}$ to be ≲1% of the ${\boldsymbol{P}}$₀ for all ${\boldsymbol{P}}$.⁹ The ratio of radiation to gas pressure ${\beta }_{r}\equiv {a}_{{\rm{R}}}{P}^{3}/(3{\rho }^{4})$, and the optical depth per wavelength $\tau \equiv \kappa \rho L$.

4.4.1. Radiation-modified Slow Mode

For a magnetized fluid in the presence of radiation, the MHD modes are damped in a similar fashion to the radiation hydrodynamic case. We first consider the radiation-modified slow mode solution. We set γ = 5/3, ρ = 1, u = 0.01, ${B}_{0}=\sqrt{5/6}$, β_r = 1, and τ = 20 and evolve the initial conditions in bhlight to t_f = 2.5, nearly half an e-folding time. The eigenmode is given in Table 1. Expected convergence at t_f in the average number of extant superphotons N_s for Monte Carlo-dominated error is shown in Figure 6.

Zoom In Zoom Out Reset image size

Table 1.  Radiation-modified Slow Mode

ω	−0.155954250795 + 0.506371984839i
δρ	0.992522043854
δu	0.0115437955397 + 0.00253571930238i
δu1	−0.0799889439467 − 0.024635280384i
δu2	−0.0804556011602 − 0.0252291311891i
δB2	−0.00672014035309 − 0.00465692766557i
δE	0.0129233747759 + 0.0201394332108i
δF1	0.00205715652365 − 0.00136719504591i
δF2	−3.27455464963 × 10−5 + 5.02957595074 × 10−5i

Download table as:  ASCIITypeset image

4.4.2. Radiation-modified Fast Mode

We now consider the radiation-modified fast mode solution, for a relativistic equilbrium. We set γ = 4/3, ρ = 1, u = 10, ${B}_{0}=\sqrt{5/6}$, β_r = 1, and τ = 20 and evolve the initial conditions in bhlight to t_f = 1.7, approximately a wave period. Note that radiation damping is not significant during this time. The eigenmode is given in Table 2. Expected convergence at t_f in the average number of extant superphotons N_s for Monte Carlo-dominated error is shown in Figure 7.

Zoom In Zoom Out Reset image size

Table 2.  Radiation-modified Fast Mode

ω	−0.000695187092855 − 3.6761842859i
δρ	0.0528655837266 − 0.000203781373769i
δu	0.704834313006 − 0.00203965222393i
δu1	0.0309307265333 − 0.000125078175647i
δu2	−0.00191512771962 + 0.000183786243559i
δB2	0.0512475714061 − 0.000472212042911i
δE	0.704838422737
δF1	−4.41864621469 × 10−5 + 0.00198608271502i
δF2	0.000397994468532 − 0.00462049532989i

Download table as:  ASCIITypeset image

Su & Olson (1996) found a solution in terms of integrals to the coupled energy balance and radiative transport equations in the diffusion approximation for a semi-infinite slab of static, initially cold fluid (with heat capacity ${c}_{v}={\mathcal{C}}{T}^{3}$) with gray absorption coefficient α, and a Marshak (isotropic incident radiation) condition at the left boundary with incident flux F. The solution is given in terms of dimensionless gas and radiation energy densities u and v, respectively, versus the dimensionless spatial coordinate $x=\sqrt{3}\alpha z$ and dimensionless time coordinate $t=4{a}_{{\rm{R}}}c\alpha t\prime /{\mathcal{C}}$, where a_R is the radiation constant. u and v are defined as:

$$\begin{eqnarray}&&u(x,t)\equiv \left(\displaystyle \frac{c}{4}\right)\left(\displaystyle \frac{E(x,t)}{F}\right),\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&v(x,t)\equiv \left(\displaystyle \frac{c}{4}\right)\left(\displaystyle \frac{{a}_{{\rm{R}}}{T}^{4}(x,t)}{F}\right),\end{eqnarray} \tag{ 58 }$$

where E is the radiation energy density.

In replicating this solution with bhlight, we adopt parameters such that the solution remains optically thick, without being so optically thick across a grid zone that our Monte Carlo transport scheme fails. We use 1024 grid zones on $x=[0,50\sqrt{3}]$ and evolve the system to t = 500 (when the solution is approximately in equilibrium). Results are shown in Figure 8.

Zoom In Zoom Out Reset image size

We now turn to dynamical tests: here, 1D radiative shocks. We will use nearly the same suite of relativistic, radiative shocks considered by Farris et al. (2008) in testing their GRRMHD code, which was based on a nonequilibrium Eddington closure. Because we solve the full transfer equation, we expect disagreement on scales of order the photon mean free path.

The Farris et al. tests assume a gray opacity κ and are set in Minkowski space. The solutions are described by ρ₀, the x component of the fluid four velocity u^x, the gas pressure P, and the comoving-frame radiation energy density E and x-component of the radiation flux four-vector F^x. The latter vanishes far from the shock. We consider only the shock-frame (unboosted) version of the tests, initialized as shock tubes. All tests are purely hydrodynamic: the magnetic field plays no role.

Our tests are identical to those given in Farris et al. (2008) except that we modify case (4), which has radiation pressure a factor of 10 larger than gas pressure upstream from the shock. Due to this large radiation pressure, bhlight cannot integrate this case stably with the numerical resources available to us. Instead, we set the upstream radiation pressure equal to the gas pressure, and call this case (4a). The shock parameters are listed in Table 3. Units are such that c = 1; a_R (equivalently, ${\rm{\hslash }}$) is determined by enforcing thermal equilibrium $E={a}_{{\rm{R}}}{(P/{\rho }_{0})}^{4}$ far from the shock.

Case 1 is a nonrelativistic strong shock with gas pressure much greater than radiation pressure, and consequently the fluid variable profiles resemble a nonrelativistic shock. bhlight output and the analytic solution are shown in Figure 9. Correspondence is good except for small deviations in the radiation variables near the shock interface, as expected for our full transfer solution.

Zoom In Zoom Out Reset image size

Case 2 is a mildly relativistic shock with somewhat larger radiation pressure than in Case 1. bhlight output and the analytic solution are shown in Figure 10. The profiles of E and F^x show qualitative differences between the bhlight result and the analytic solution: a discontinuity in the analytic solution does not appear for the case of full transfer. This is unsurprising given the approximate nature of the Farris solutions. Note that for this case the approximate Eddington solution contains an unphysical discontinuity in the coordinate frame radiation energy density.

Zoom In Zoom Out Reset image size

Case 3 is a highly relativistic shock with dynamically important radiation field. bhlight output and the analytic solution are shown in Figure 11. We find significant differences within a few photon mean free paths of the shock, particularly for the radiation flux. Figure 12 shows the expected N^−1/2 self-convergence in the radiation variables. A similar self-convergence trend also appears in the fluid variables until grid resolution becomes the dominant source of error.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Case 4a is a modestly relativistic wave with upstream radiation and gas pressure nearly equal. bhlight output and the analytic solution are shown in Figure 13. We find good agreement in all variables, despite the relatively strong radiation field. Note also that even with a large number of samples it is difficult to suppress noise in F^x when it is much smaller than E.

Zoom In Zoom Out Reset image size

Next we turn to a general relativistic equilibrium test. Consider a Schwarzschild black hole surrounded by a static atmosphere. The atmosphere is bounded by static, concentric spherical shells at ${r}_{i}\gt 2{GM}/{c}^{2}$ and ${r}_{o}\gt {r}_{i}$, and is in radiative equilibrium. The shells are reflecting boundaries and exchange no heat with the atmosphere, which has a gray opacity κ and adiabatic index γ = 5/3.

In the Newtonian limit radiative conduction would drive the atmosphere toward T = const. In a relativistic atmosphere the gravitational field causes the atmosphere to come into a different equilibrium in which redshifted temperature, ${T}_{\infty }\equiv T\sqrt{-{g}_{00}}$, is constant.

In the Eddington approximation the atmospheric structure is determined by the conditions ${T}_{\infty }=\mathrm{const}$, and mechanical equilibrium, ${T}_{;r}^{{rr}}=0$. Note that any radiation source terms vanish in thermal equilibrium; in fact, the Eddington approximation solution turns out to be an exact solution to the transfer equation, as we show in Appendix B.

Once the inner and outer radii are specified, there are three dimensionless parameters that describe the solution (although their interpretation is purely Newtonian): the ratio of the atmospheric scale height at the inner boundary to the local radius $h\equiv \gamma {k}_{{\rm{B}}}{T}_{i}{r}_{i}/(\mu {m}_{p}{GM})$ (where μ is the mean molecular weight), the ratio of radiation pressure to gas pressure at the inner boundary ${\beta }_{r}\equiv \mu {m}_{p}{a}_{{\rm{R}}}{T}_{i}^{3}/(3{\rho }_{i}{k}_{{\rm{B}}})$, and the characteristic optical depth $\tau \equiv \kappa {\rho }_{i}({r}_{o}-{r}_{i})$.

We set ${r}_{i}=3{GM}/{c}^{2}$, ${r}_{o}=20{GM}/{c}^{2}$. We set h ≈ 2.66, β_r ≈ 0.23, and τ ≈ 5.0. We use a 1D domain with 128 zones; the errors are dominated by Monte Carlo noise. The solution is shown in Figure 14. We also find the comoving spectrum at fiducial radius r ≈ r_i to match our expectation for thermalized radiation (with free vertical scale), as shown in Figure 15. Evidently the redshifted temperature is indeed constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we report a preliminary investigation of spherically symmetric accretion onto a Schwarzschild black hole with radiative coupling. The fully dynamical problem has been studied analytically (although not with a full transport solution, frequency dependence, or magnetic fields) by Vitello (1984), Park (1990), and Nobili et al. (1991). Frequency-integrated numerical studies have also been performed in Fragile et al. (2012), Roedig et al. (2012), and McKinney et al. (2014) over many GM/c². Here we obtain a full transport solution with frequency dependence and magnetic fields. We include synchrotron emission, synchrotron absorption, and Compton scattering, with associated heating and cooling.

We set M = 6.6 × 10⁹ M_⊙ and describe the spacetime with MKS coordinates (Section 3.5). All simulations are performed in 1D, with $r\in [1.5{\text{}}M,50{\text{}}M]$ for 64 zones. As initial conditions for the fluid, we adopt the nonradiative Bondi solution of Hawley et al. (1984; which we hold constant at the outer boundary), except we set γ = 13/9 and place the sonic point at r = 200M. We vary the accretion rate by varying the density (or equivalently the mass unit ${\mathcal{M}}$) of the flow.

The magnetic field is initialized as a radial, monopolar field, which has no effect on the fluid motion. We set ${B}^{1}=\alpha /{r}^{3}$, where α is chosen such that $\beta \equiv 2\;P/{b}^{2}\approx 130$ at r = 2M in the initial (GRMHD) conditions. No radiation is present initially. The radiation is allowed to flow freely out of the computational domain at the radial boundaries, with no inflow of radiation.

The luminosity L(r), evaluated once the flow has settled and become almost time-independent, is

$$\begin{eqnarray}&&L(r)\equiv \int \sqrt{-g}\;{{dx}}^{2}{{dx}}^{3}{R}_{0}^{1}.\end{eqnarray} \tag{ 59 }$$

We average $L(r)$ over radial shells between $r=30{GM}/{c}^{2}$ and $r=50{GM}/{c}^{2}$ to obtain an average luminosity L. This is equivalent to time-averaging at a single radius.

We perform five calculations at different accretion rates. We evolve each system until the luminosity becomes stable. At high accretion rate this can take as long as $400{GM}/{c}^{3}$. The characteristic number density for each simulation, mass accretion rate

$$\begin{eqnarray}&&\dot{M}\equiv -\int \sqrt{-g}\;{{dx}}^{2}{{dx}}^{3}\rho {u}^{1}\end{eqnarray} \tag{ 60 }$$

(evaluated at the event horizon), and L are given in Table 4. The resulting profiles of these five cases, along with that of Case 5 with Compton scattering disabled, are shown in Figure 16. The lack of cooling in the pure synchrotron case indicates that Compton cooling dominates over synchrotron cooling for the highest accretion rate model. The relationship between luminosity and accretion rate when Compton scattering is active is shown in Figure 17.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Parameters for Radiating Bondi Accretion

Case	n	$\dot{m}$	L
	(cm−3)		LEdd
1	3.0 × 106	4.01 × 10−7	2.03 × 10−14
2	3.0 × 107	4.01 × 10−6	1.46 × 10−12
3	3.0 × 108	4.01 × 10−5	2.08 × 10−10
4	3.0 × 109	4.01 × 10−4	5.81 × 10−7
5	3.0 × 1010	4.01 × 10−3	3.05 × 10−4

Download table as:  ASCIITypeset image

We have checked self-convergence of a solution with n = 1.51 × 10¹⁰ cm⁻³ (between Cases 4 and 5 in Table 4) in steady state at t = 200M. The expected convergence behavior for Monte Carlo-dominated error is shown in Figure 18.

Zoom In Zoom Out Reset image size

Our models show that Compton cooling is important for Bondi accretion near the Eddington rate, and that—for our assumed field configuration—synchrotron cooling is comparatively unimportant. Although here we find appreciable cooling only close to the Eddington rate, we expect that for near-Keplerian accretion flows Compton scattering will become dynamically important at lower accretion rates, as individual fluid elements will have more time to cool (their radial velocity is lower) before accreting onto the black hole.

As another preliminary application of bhlight we consider the effect of radiation on an intermediate accretion rate black hole accretion flow. Recall that for systems with L ∼ 10⁻⁹ L_Edd (like Sgr A*) radiation will have little effect since the cooling timescale is long compared to the accretion timescale. This is the classical RIAF regime. For systems with L ∼ 10⁻⁵ L_Edd (like M87) the effect of radiation depends on temperature and the distribution function of the electrons, but in certain circumstances radiation interactions—especially Compton cooling—can cool the flow on timescales comparable to or shorter than the accretion timescale.

Here we consider an axisymmetric accretion flow with parameters inspired by M87. We set black hole mass M = 6.6 × 10⁹ M_⊙ following Gebhardt et al. (2011) and dimensionless spin a_* = 0.9375, and we adjust ${\mathcal{M}}$ so that $\dot{m}\approx 6.3\times {10}^{-6}$ (Kuo et al. 2014). Although the distribution functions for ions and electrons in M87-like systems probably contain multiple components, we adopt thermal distributions for both and set the ion and electron temperatures T_i = T_e. T_i/T_e may strongly affect the dynamics of this system; essentially, it controls the cooling rate. We set the adiabatic index γ = 13/9, appropriate for ionized hydrogen when ${{kT}}_{i}\ll {m}_{p}{c}^{2}$ and ${{kT}}_{{\rm{e}}}\gg {m}_{{\rm{e}}}{c}^{2}$. We include synchrotron emission for a relativistic, thermal distribution of electrons (emissivity given by Equation (72) of Leung et al. 2011), thermal absorption, and Compton scattering. Bremsstrahlung is neglected, as it is a small correction to the emissivity within ∼10³ M of the black hole for such intermediate accretion rate systems (e.g., Narayan & Yi 1995).

The initial conditions are a Fishbone–Moncrief torus (Fishbone & Moncrief 1976) with an inner radius at $6{GM}/{c}^{2}$ and pressure maximum at $12{GM}/{c}^{2}$. We extend this configuration by adding a weak poloidal magnetic field that follows isodensity contours, using the same procedure as in McKinney & Gammie (2004), but with the vector potential modified by a $\mathrm{cos}\theta$ factor to produce a two-loop configuration. The field strength is normalized so that β has a global minimum value of 100. Small perturbations are applied to the internal energy to efficiently initiate the MRI. No radiation is present in the initial conditions.

We use the MKS coordinates (see Section 3.5), with r₀ = 0 and h_s = 0.3. The grid runs from $r=0.98(1+\sqrt{1-{a}_{*}^{2}}){GM}/{c}^{2}$ to $r=40{GM}/{c}^{2}$, and from θ = 0 to θ = π. We impose outflow boundary conditions on both the fluid and radiation at the inner and outer radial boundaries, and reflecting polar boundary conditions. In this instance we adopt the piecewise parabolic method for reconstruction. We set the CFL number to 0.7, and evolve the system until t = 2000M.

To accurately sample cooling due to scattering events, we bias the scattering via the method given in Section 3.2.7 such that ${b}_{{\rm{s}}}\propto {{\rm{\Theta }}}_{{\rm{e}}}^{2}$. To determine the absolute number of superphotons in steady state, we require only that the bolometric light curve be satisfactorily resolved (here, ∼1.1 × 10⁶ superphotons at any one time). In the future we will more carefully consider the resolution requirement for these models.

We find qualitative differences between the radiative torus simulated in bhlight and the same model evolved with ideal GRMHD due mainly to Compton cooling in the hot, dense regions of the flow. Figure 19 shows fluid temperature and comoving radiation energy density at $t=2000M$. Figure 20 shows density contours at t = 2000M for the two models. Figure 21 shows shell-averaged density-weighted temperature,

$$\begin{eqnarray}&lt;{{\rm{\Theta }}}_{{\rm{e}}}&gt;\equiv \displaystyle \frac{\int \sqrt{-g}{{dx}}^{2}{{dx}}^{3}\rho {{\rm{\Theta }}}_{{\rm{e}}}}{\int \sqrt{-g}{{dx}}^{2}{{dx}}^{3}\rho },\end{eqnarray} \tag{ 61 }$$

as a function of radius for both bhlight and ideal GRMHD calculations; evidently the plasma is significantly cooler in the model with radiation interactions. Notice, however, that our single temperature model may be having an unrealistically strong dynamical effect; at these accretion rates the ions are likely imperfectly coupled to the electrons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We evaluate the luminosity L, Equation (59), at large radius by integrating over all zones in a spherical shell. Figure 22 shows L and the radiative efficiency $\eta =L/\dot{M}{c}^{2}$. The mean η we find between t = 1750M and t = 2000M, $\langle \eta \rangle =0.51$, where

$$\begin{eqnarray}&&\langle \eta \rangle \equiv \displaystyle \frac{\int {dt}\;L}{\int {dt}\;\dot{M}{c}^{2}},\end{eqnarray} \tag{ 62 }$$

is high compared to the thin disk value, $\approx 0.18$ (for this a_*). Most of this energy is extracted by Compton scattering at $r\sim 10-15{GM}/{c}^{2}$. This high efficiency is a consequence of the flow having not yet reached steady state.

Zoom In Zoom Out Reset image size

We plan to explore radiative flows at intermediate accretion rates with more sophisticated treatments of the electron physics in future work.

The static atmosphere test in Section 4.7 uses reflecting boundaries that are aligned with the coordinates. How should one apply the reflecting boundary conditions to the radiation field? This is not trivial in Kerr–Schild coordinates where the naive approach of simply changing the sign of the radial component of the wavevector is wrong, because the radial component of the shift does not vanish.

When a photon crosses the boundary, we build an orthonormal tetrad with time component ${e}_{(t)}^{\mu }=({u}^{t},0,0,0)$ (i.e., the tetrad is stationary in the coordinate, and hence the boundary, frame) and one spatial component that is normal to the boundary. We transform the wavevector to the tetrad frame, reverse the normal component of the wavevector, and transform back to the coordinate frame.

For problems with fluid inflow across the boundary (e.g., the relativistic shocks in Section 4.6) the fluid advects a thermal radiation field with it across the boundary. How should one sample the incoming photons on the boundary? The problem is that one is sampling a flux rather than the distribution function itself.

We have found that the simplest procedure is to sample the distribution function and multiply the weights in the sample by $\mathrm{cos}\theta$, where θ is the angle between the wavevector and the normal to the boundary.