grmonty: A MONTE CARLO CODE FOR RELATIVISTIC RADIATIVE TRANSPORT

There is wide interest in calculating the emergent radiation from relativistic astrophysical sources, including accreting black holes, accreting neutron stars, and relativistic blast waves. A variety of methods for solving the radiative transfer problem in these sources have been developed over the last few decades (e.g., Pozdynakov et al. 1983; Górecki & Wilczewski 1984; Hauschildt & Wehrse 1991; Carrigan & Katz 1992; Coppi et al. 1993; Stern et al. 1995; Poutanen & Svensson 1996; Zane et al. 1996; Dove et al. 1997; Böttcher & Liang 2001; Schnittman et al. 2006; Noble et al. 2007; Wu et al. 2008), some based on Monte Carlo schemes. Few of these schemes take full account of relativistic effects in the source, however, and this is crucial in estimating the spectra of hot plasma deep in a gravitational potential well, or highly relativistic blast waves. Monte Carlo transport of radiation in accretion flows around compact objects has been considered by (e.g., Schnittman & Krolik 2009; Schnittman 2006; Yao et al. 2005; Böttcher et al. 2003; Böttcher & Liang 2001; Laurent & Titarchuk 1999; Agol & Blaes 1996; Stern et al. 1995). Among others, Cullen (2001), Molnar & Birkinshaw (1999), Hua (1997), Górecki & Wilczewski (1984), and Pozdynakov et al. (1983) give more general discussions of Monte Carlo radiative transfer techniques.

We were motivated by efforts to model the radio source and black hole candidate Sgr A*. Our interest in this source drove us to develop a numerical scheme that could accurately calculate spectra of a relativistic source in which the plasma properties (velocity, density, magnetic field strength, and temperature) were specified by a separate model—that is, sources in which radiation plays a negligible role in the dynamics and energetics. The result, a Monte Carlo scheme called grmonty, is described in this paper. The spirit of our calculation is to obtain an accurate spectrum with as few approximations as possible. To this end, we treat Compton scattering with no expansions in v/c, and allow for general angle-dependent emission and absorption (we specialize to thermal synchrotron in this work).

In designing grmonty, our philosophy has been to maximize the physical transparency and minimize the length of the code, occasionally at the cost of reduced performance. Sometimes simplicity and efficiency are in harmony. We chose to directly integrate the geodesic equation rather than using a scheme that relies on integrability of geodesics in the Kerr metric. We will show that for radiative transfer problems where many points are required along each geodesic, direct integration is not only simpler and easier to modify, but also faster.

Our paper is organized as follows. In Section 2, we describe how we sample emission, and in Section 3 we describe how we track photons along geodesics. Evolution of superphoton weights under absorption is described in Section 4, and sampling of scattered photons is discussed in Section 5. Photons at large distance from the source must be sampled and assembled into spectra; this is described in Section 6. The code has been extensively verified; we describe the tests in Section 7. Section 8 describes a sample calculation and Section 9 summarizes our results.

Throughout this paper we assume that there is an underlying model that can be queried to supply the rest-mass density ρ₀, the internal energy u, the four-velocity u^μ, and magnetic field four-vector b^μ for the radiating plasma. Usually we expect the model to be supplied by a numerical simulation in a coordinate basis x^μ.

Emission in grmonty is treated by sampling the emitted photon field. The samples, here called "superphotons" (also "photon packets"), have weight w, coordinates x^μ, and wave vector k^μ. The weight w ≫ 1 is a pure number that represents the ratio of photons to superphotons: dN = wdN_s (N_s≡ number of superphotons, N≡ number of photons). In our models the weight is a function of the emitting plasma frame frequency ν and nothing else. The coordinates are typically in model units (e.g., for a black hole accretion flow calculation, length unit L = GM/c²), and the components of k^μ are given in units of m_ec².

How should superphotons be distributed over x^μ and k^μ? The initial superphoton momentum can be described in an orthonormal tetrad basis e^μ_(a) that is attached to the plasma, so that e^μ_(t) = u^μ (μ is the coordinate index, and (a) is the index associated with the tetrad basis, raised and lowered using the Minkowski metric). In the tetrad basis k^(a) is specified by frequency ν and spatial direction unit vector $\hat{\bf n}$ that is contained within the solid angle dΩ. The probability distribution for superphotons is then

$$\begin{equation} \frac{1}{{\sqrt{-g}}}\, \frac{d N_s}{d^3x dt \, d\nu d\Omega } = \frac{1}{w{\sqrt{-g}}} \, \frac{d N}{d^3x dt \, d\nu d\Omega } = \frac{1}{w} \frac{j_\nu }{h\nu }, \end{equation} \tag{ 1 }$$

where j_ν is the emissivity (always defined in the plasma frame), since ${\sqrt{-g}}\, d^3x dt$ is invariant (meaning coordinate invariant). In a time interval Δt, we expect to create

$$\begin{equation} N_{s,{\rm tot}} = \Delta t \int {\sqrt{-g}}\; d^3x \; d\nu \; d\Omega \frac{1}{w}\frac{j_\nu }{h\nu } \end{equation} \tag{ 2 }$$

superphotons over the entire model volume. The total computational effort is proportional to N_s,tot, so we control the computational effort by scaling the weights.

How should we distribute superphotons over the volume? grmonty subdivides the model volume into volume elements ("zones") of size Δ³x (e.g., in Boyer-Lindquist t, r, θ, ϕ coordinates for the Kerr metric, ΔrΔθΔϕ). For zone i we expect

$$\begin{equation} N_{s,i} = \Delta t \Delta ^3x {\sqrt{-g}}\int \! d\nu \; d\Omega \; \frac{1}{w}\, \frac{j_\nu }{h\nu } \end{equation} \tag{ 3 }$$

to be created in time Δt. We create N_s,i superphotons at the center of zone i. Fractional N_s,i are dealt with by rejection sampling.

The momentum-space (wave vector) piece of the probability distribution (1) can be sampled by techniques outlined below to give ν and $\hat{\bf n}$. With x^μ, ν and $\hat{\bf n}$ in hand, we can construct k^(a) and, finally, transform to the coordinate basis e^μ_(a)k^(a) = k^μ.

The remaining ingredients in the sampling procedure are the emissivity, the orthonormal tetrads, and the sampling procedures for ν and $\hat{\bf n}$.

2.1. Emissivity

grmonty depends on the emissivity only through functions that specify j_ν and ∫dνdΩj_ν/ν, so it is straightforward to include any emission/absorption process.

In our target problem, the only source of superphotons is thermal synchrotron emission at dimensionless temperature Θ_e ≡ kT_e/(m_ec²). P. K. Leung & C. F. Gammie (2009, in preparation) show that, for Θ_e ≳ 0.5,

$$\begin{equation} j_\nu (\nu,\theta) \simeq \frac{\sqrt{2}\pi e^2 n_e\nu _s}{3 c K_2(\Theta _e^{-1})} (X^{1/2} + 2^{11/12} X^{1/6})^2 \exp (-X^{1/3})\quad \end{equation} \tag{ 4a }$$

$$\begin{equation} X \equiv \frac{\nu }{\nu _s} \end{equation} \tag{ 4b }$$

$$\begin{equation} \nu _s \equiv \frac{2}{9} \left(\frac{e B}{2\pi m_e c}\right) \Theta _e^2 \sin \theta \end{equation} \tag{ 4c }$$

where K₂ is the modified Bessel function of the second kind, n_e is the number density of electrons, B is the magnetic field strength, and θ is the angle between the wave vector and magnetic field. For large Θ_e, K₂(Θ⁻¹_e) ≃ 2Θ²_e, but for Θ_e ≲ 1 better agreement with the emissivities of P. K. Leung & C. F. Gammie (2009, in preparation) is obtained if K₂ is evaluated directly. Since the emissivity must be evaluated many times, it is most efficient to precompute K₂(Θ⁻¹_e) at the beginning of the calculation and store the results in a table.

2.2. Orthonormal Tetrads

The wave vector sampling is done in an orthonormal tetrad attached to the fluid. We construct the orthonormal tetrad e^μ_(a) using numerical Gram–Schmidt orthogonalization. Here μ is the coordinate index, and (a) is the index associated with the tetrad basis, which is raised and lowered using the Minkowski metric.

We set e^μ₍₀₎ = u^μ (u^μ≡ plasma four-velocity), and then use b^μ, the magnetic field four-vector, as the first trial vector (this is numerically convenient since we will want to orient wave vectors with respect to the magnetic field; if b^μ = 0 then we use a default, radius-aligned, trial vector). Thus e^μ₍₁₎ = NORM(b^μ − e^μ₍₀₎(e^ν₍₀₎b_ν)) (NORM: normalize). The process is repeated with additional trial vectors to create a full tetrad basis.

The tetrad-to-coordinate basis transformation is

$$\begin{equation} k^\mu = e^\mu _{(a)} k^{(a)} \end{equation} \tag{ 5 }$$

and the coordinate to tetrad transformation is

$$\begin{equation} k^{(a)} = e^{(a)}_\mu k^\mu. \end{equation} \tag{ 6 }$$

With these transformations in hand, we can construct the superphoton wave vector in the orthonormal frame and then transform it to the coordinate basis.

2.3. Wave Vector Sampling Procedure

2.3.1. Photon Energy

Within zone i superphotons are distributed over frequency according to the distribution

$$\begin{equation} \frac{dN_{s,i}}{d\ln \nu } = \Delta t \Delta ^3 x {\sqrt{-g}}\; \frac{1}{h\, w} \; \int d\Omega \; j_\nu. \end{equation} \tag{ 7 }$$

We sample the distribution only between a minimum and maximum frequency ν_min and ν_max. These must be chosen so that no significant emission is omitted from the final spectrum.

We distribute superphotons over frequency by rejection sampling. For simplicity, we use a constant envelope function equal to the maximum of Equation (7) (for zone i). Thus we draw tentative values uniformly in ln ν from ν_min to ν_max. The efficiency of the sampling procedure is the ratio of the areas under the distribution and envelope, so if the distribution given in Equation (7) is sharply peaked this technique can be inefficient.

In practice, we choose a tentative frequency ν₀ = exp(r₁ln ν_max/ν_min + ln ν_min), where r₁ is drawn from a uniform distribution on [0,1) (we use the Mersenne twister random number generator from the GNU Scientific Library, hereafter GSL). A second number r₂ is drawn from [0,1) and the process is repeated until

$$\begin{equation} r_2 < {\left.\frac{dN_{s,i}}{d\ln \nu }\right|}_{\nu _0}\bigg/ {\rm MAX}\left(\frac{dN_{s,i}}{d\ln \nu }\right). \end{equation} \tag{ 8 }$$

The efficiency of this process is ∼15%, but the cost is small compared to the total cost of grmonty.

2.3.2. Photon Direction

The superphoton direction $\hat{\bf n}$ is described by polar coordinates θ and ϕ in the tetrad frame, where θ is the angle between the spatial part of the wave vector and the magnetic field. The colatitude θ is obtained by rejection sampling: a tentative value for θ is obtained by drawing μ = cos θ from a uniform distribution on [ − 1, 1), a second number r is drawn from a uniform distribution on [0, 1), and θ is accepted if

$$\begin{equation} r < \frac{j_\nu (\theta)}{j_\nu (\pi /2)} \end{equation} \tag{ 9 }$$

(this procedure is specific to the synchrotron emissivity). The efficiency of this scheme is problem dependent; for our target application the efficiency is ∼65%. Finally, ϕ is drawn from a uniform distribution on [0, 2π).

2.3.3. Transformation to Coordinate Frame

Once θ, ϕ, and = hν/m_ec² are selected, the wave vector is completely specified in the orthonormal tetrad frame:

$$\begin{equation} k^{(0)} = {\epsilon } \end{equation} \tag{ 10a }$$

$$\begin{equation} k^{(1)} = {\epsilon }\cos \theta \end{equation} \tag{ 10b }$$

$$\begin{equation} k^{(2)} = {\epsilon }\sin \theta \cos \phi \end{equation} \tag{ 10c }$$

$$\begin{equation} k^{(3)} = {\epsilon }\sin \theta \sin \phi, \end{equation} \tag{ 10d }$$

and the wave vector in the coordinate frame is k^μ = e^μ_(a)k^(a).

General relativistic radiative transfer differs from conventional radiative transfer in Minkowski space in that photon trajectories are no longer trivial; photons move along geodesics. Tracking geodesics is a significant computational expense in grmonty.

The governing equations for a photon trajectory are

$$\begin{equation} \frac{d x^\alpha }{d\lambda } = k^\alpha \end{equation} \tag{ 11 }$$

which defines λ, the affine parameter, the geodesic equation

$$\begin{equation} \frac{d k^\alpha }{d\lambda } = -\Gamma ^\alpha _{\mu \nu } k^\mu k^\nu, \end{equation} \tag{ 12 }$$

and the definition of the connection coefficients

$$\begin{equation} \Gamma ^\alpha _{\mu \nu } = \frac{1}{2} g^{\alpha \gamma } (g_{\gamma \mu,\nu } + g_{\gamma \nu,\mu } - g_{\mu \nu,\gamma }) \end{equation} \tag{ 13 }$$

in a coordinate basis.

We assume nothing about the metric, so it is easy to change coordinate systems and even to extend the code to dynamical spacetimes. Nevertheless, our main application—to black hole accretion flows—is in the Kerr metric, where geodesics are integrable. The four constants of the motion are the energy-at-infinity E (in Boyer–Lindquist coordinates t, r, θ, ϕ, E = −k_t), the angular momentum l = k_ϕ, Carter's constant $Q = k_\theta ^2 + k_\phi ^2 \cot ^2\theta - a^2 k_t^2 \cos ^2\theta$ (see Carter 1968), and the condition that k^α be null: k^μk_μ = 0, equivalent to the dispersion relation for photons in vacuo: ω² = c²k². These four constants of the motion can be used to quasi-analytically obtain x^μ and k^μ in terms of an initial (or final) position and wave vector (see, e.g., Rauch & Blandford 1994; Beckwith & Done 2005; Dexter & Agol 2009).

The integrability of geodesics in the Kerr metric would appear to provide an opportunity for significant computational economies. We show below, however, that direct integration of Equations (11) and (12) is not only simpler and more flexible but also faster than at least one implementation of an integral-based technique.

Which ODE integration algorithm is best for the geodesic equation? If only a few coordinate evaluations are required over the entire geodesic then a high-order scheme is optimal. For example, we have found that the embedded Runge–Kutta Prince–Dorman method available in GSL is fast and accurate; it can easily be made to conserve the integrals of motion to machine precision. Many coordinate evaluations are required, however, when integrating the equation of radiative transfer, as grmonty does, along superphoton trajectories. A second-order scheme can then provide the required accuracy at minimal cost.

Evaluating the connection coefficients is expensive, so we want to choose a scheme that minimizes the number of evaluations. The velocity Verlet algorithm, which for the geodesic equation is

$$\begin{equation} x^{\alpha }_{n+1}=x^{\alpha }_n+k^{\alpha }_n\Delta \lambda +{{1}\over {2}}\left({{dk^{\alpha }}\over {d\lambda }}\right)_n(\Delta \lambda)^2 \end{equation} \tag{ 14a }$$

$$\begin{equation} k^{\alpha }_{n+1,p}=k^{\alpha }_n+\left({{dk^{\alpha }}\over {d\lambda }}\right)_n\Delta \lambda \end{equation} \tag{ 14b }$$

$$\begin{equation} \left({{dk^{\alpha }}\over {d\lambda }}\right)_{n+1}=-\Gamma ^{\alpha }_{\mu \nu }(\mathbf {x}_{n+1})k^{\mu }_{n+1,p}k^{\nu }_{n+1,p} \end{equation} \tag{ 14c }$$

$$\begin{equation} k^{\alpha }_{n+1}=k^{\alpha }_n + {{1}\over {2}}\left(\left({{dk^{\alpha }}\over {d\lambda }}\right)_n + \left({{dk^{\alpha }}\over {d\lambda }}\right)_{n+1}\right)(\Delta \lambda), \end{equation} \tag{ 14d }$$

requires only one evaluation of the connection coefficients per step. Accuracy can be improved by using the result of Equation (14d) to recompute the derivative Equation (14c) with k^μ_n+1,p = k^μ_n+1 and then reevaluating Equation (14d). This process is repeated until the change in the wave vector between estimates is below some tolerance. In grmonty , we continue this iteration until the fractional change is less than 10⁻³ (typically only once or twice). This does not require any additional evaluations of Γ^α_μν. Very rarely we find that this iteration fails to converge, and then grmonty defaults to taking the step with a classical fourth-order Runge–Kutta technique.

How fast and accurate is our geodesic integration scheme? We propose the following benchmark. Consider a point on a direct, circular, marginally stable orbit in the equatorial plane of a black hole with spin a/M = 1 − 2⁻⁴ = 0.9375 that emits radiation isotropically in its rest frame. Sample the emitted photons (in a Monte Carlo sense; the analytic circular orbit orthonormal tetrads available in Bardeen et al. 1972 are useful for constructing the initial wave vectors) and track them until they cross the horizon or reach rc²/(GM) = 100 (r is the Boyer–Lindquist or Kerr–Schild radial coordinate). Figure 1 shows as dots a representative sample of photon geodesics from grmonty in the coordinate frame, illustrating the effects of relativistic beaming, lensing, and frame dragging.

Zoom In Zoom Out Reset image size

Second-order convergence of the velocity Verlet integration scheme is demonstrated in Figure 2, which plots the average fractional error in E, l, and Q as a function of a step-size parameter ε. We typically set ε = 0.04 as a compromise between performance and accuracy; the average fractional errors at the end of the integrations are ∼2 × 10⁻³, ∼4 × 10⁻², and ∼8 × 10⁻² for E, l, and Q, respectively. We have verified that this choice makes geodesic tracking errors subdominant in the error budget for the overall spectrum.

Zoom In Zoom Out Reset image size

With ε = 0.04, grmonty integrates ∼16, 700geodesicss⁻¹ on a single core of an Intel Xeon model E5430. If we use fourth-order Runge–Kutta exclusively so that the error in E, l, and Q is ∼1000 times smaller, then the speed is ∼6200geodesicss⁻¹. If we use the Runge–Kutta Prince–Dorman method in GSL with ε = 0.04 the fraction error is ∼10⁻¹⁰ and the speed is ∼1700geodesicss⁻¹. These results can be compared to the publicly available integral-based geokerr code of Dexter & Agol (2009), whose geodesics are shown as the (more accurate) solid lines in Figure 1. If we use geokerr to sample each geodesic the same number of times as grmonty (∼180), then on the same machine geokerr runs at ∼1000geodesicss⁻¹. It is possible that other implementations of an integral-of-motion-based geodesic tracker could be faster.

If only the initial and final states of the photon are required, we find that geokerr computes ∼77, 000geodesicss⁻¹. The adaptive Runge–Kutta Cash–Karp integrator in GSL computes ∼34, 500geodesicss⁻¹ with fractional error ∼10⁻³.

grmonty treats absorption deterministically. We begin with the radiative transfer equation written in the covariant form

$$\begin{equation} \frac{1}{{\mathcal C}} \frac{d}{d\lambda } \left(\frac{I_\nu }{\nu ^3}\right) = \left(\frac{j_\nu }{\nu ^2}\right) - (\nu \alpha _{\nu,a}) \left(\frac{I_\nu }{\nu ^3}\right). \end{equation} \tag{ 15 }$$

(see Mihalas & Mihalas 1984). Here I_ν is specific intensity and α_ν,a is the absorption coefficient (which is always evaluated in the fluid frame). The absorption coefficient must be computed by a separate subroutine; for thermal synchrotron emission we set α_ν,a = j_ν/B_ν. ${\mathcal C}$ is a constant that depends on the units of k^μ (in grmonty, electron rest mass), ν (Hz), and the length unit L for the simulation in cgs units. For grmonty,

$$\begin{equation} {\mathcal C}= \frac{L h}{m_e c^2}. \end{equation} \tag{ 16 }$$

Each quantity in parentheses in Equation (15) is invariant; I_ν/ν³, for example, is proportional to the photon phase space density.

Since I_ν/ν³ is proportional to the number of photons moving along each ray, I_ν/ν³ ∝ w, and Equation (15) implies (ignoring emission)

$$\begin{equation} \frac{d w}{d\tau _{a}} = -w, \end{equation} \tag{ 17 }$$

where

$$\begin{equation} d\tau _{a} = (\nu \alpha _{\nu,a})\; {\mathcal C}d\lambda \end{equation} \tag{ 18 }$$

is the differential optical depth to absorption and the quantity in parentheses is the "invariant opacity." This equation we integrate with second-order accuracy

$$\begin{equation} \tau _{a}=\frac{1}{2}((\nu \alpha _{\nu,a})_n + (\nu \alpha _{\nu,a})_{n+1}) \; {\mathcal C}\Delta \lambda, \end{equation} \tag{ 19 }$$

and then set

$$\begin{equation} w_{n+1} = w_{n} e^{-\tau _{a}}. \end{equation} \tag{ 20 }$$

Since the components of k^μ are expressed in units of the electron rest-mass energy, ν = −k^μu_μm_ec²/h. Storing the invariant opacity at the end of each step saves computations since it can be reused as the beginning of the following step.

Our treatment of scattering consists of two parts: the first determines where a superphoton should scatter and the second determines the energy and direction of the scattered superphoton.

5.1. Selection of Scattering Optical Depth

When a superphoton is created or scattered, grmonty selects the scattering optical depth τ_s at which the next scattering event will take place. Scattering follows the cumulative probability distribution

$$\begin{equation} p = 1 - e^{-\tau _{s}} = \tau _{s} + O\big(\tau _{s}^2\big), \end{equation} \tag{ 21 }$$

so superphotons will experience on average τ_s scattering events when τ_s ≲ 1. In optically thin sources this would result in poor signal to noise in portions of the spectrum dominated by scattered light. To overcome this, we use the biased probability distribution

$$\begin{equation} p = 1- e^{-b\tau _s}, \end{equation} \tag{ 22 }$$

where b is a bias parameter. This technique was originally proposed by Kahn (1950) in the context of deep penetration of neutrons in radiation shielding and has since been extensively explored in the nuclear engineering literature (often referred to as exponential biasing or exponential transform). Whereas in deep penetration problems b ⩽ 1 in order to allow for sampling of radiation at high optical depths, here we set b ⩾ 1 in order to better sample scattered photons at low optical depths. Superphotons now experience on average bτ_s scattering events. Two superphotons emerge from a scattering event: the incident superphoton of weight w and a new scattered superphoton. For conservative scattering the incident superphoton has its weight reset to w(1 − 1/b) and the new superphoton has weight w/b, so that weight (photon number) is conserved.

What should we choose for the bias parameter b? The goal is to set b such that scattering is more likely to occur in regions which contribute most to the spectrum. This is an example of the more general technique of importance sampling. Typically we set the bias parameter b = MAX(1, αΘ²_e/τ_s,max) (α is a scaling factor we set to 1/〈Θ_e〉² where 〈Θ_e〉 is the volume averaged dimensionless temperature and τ_s,max is an estimated maximum scattering optical depth) to improve sampling on the high-energy side of each scattering order, which is populated by photons scattered from a high-temperature plasma. If the bias factor is too large a "chain reaction" results in an exponentially growing number of superphotons; 1/τ_s,max is an estimate of the critical bias factor in a Θ_e = 1 plasma.

We evaluate the scattering optical depth along geodesics in a manner analogous to the absorption optical depth;

$$\begin{equation} \tau _{s}=\frac{1}{2}((\nu \alpha _{\nu,s})_n + (\nu \alpha _{\nu,s})_{n+1}) \; {\mathcal C}\Delta \lambda \end{equation} \tag{ 23 }$$

is the scattering depth along a step.

5.1.1. Covariant Evaluation of Extinction Coefficient

In our applications electron scattering dominates. The general, invariant expression for the rate of binary interactions dN_ab between a population of particles dN_a and dN_b is

$$\begin{eqnarray} \frac{1}{{\sqrt{-g}}}\,\frac{dN_{ab}}{d^3x dt} &=& \frac{1}{1 + \delta _{ab}}\; \int \; \frac{d^3p_a}{{\sqrt{-g}}\,p_a^t} \; \frac{d^3p_b}{{\sqrt{-g}}\,p_b^t}\nonumber\\ &\times &\frac{dN_a}{d^3x \, d^3p_a} \; \frac{dN_b}{d^3x \, d^3p_b} \big({-}p_{a\mu } p_b^\mu\big) \; \sigma v_{ab},\quad \end{eqnarray} \tag{ 24 }$$

where δ_ab prevents double counting if a = b, d³p = dp₁dp₂dp₃, σ is the invariant cross section, and v_ab = c(1 + m²_am²_b/(−p_aμp^μ_b)²)^1/2. This is the manifestly covariant generalization of Equation (12.7) of Landau & Lifshitz (1975).

We want to use this expression to find the cross section for a photon with wave vector k^μ₀ interacting with a population of particles of mass m>0. We therefore set dN_γ/d³xd³p = δ(k^μ − k^μ₀) and, dropping the subscripts on k and p, the integral reduces to

$$\begin{equation} \frac{1}{{\sqrt{-g}}}\,\frac{dN_{m\gamma }}{d^3x dt} = \int \; \frac{d^3p}{{\sqrt{-g}}\,p^t} \; \frac{dn_m}{d^3p} \; \frac{(-k_{\mu } p^\mu)}{k^t} \; \sigma c, \end{equation} \tag{ 25 }$$

where dn_m = dN_m/d³x. In Minkowski coordinates (${\sqrt{-g}}= 1$), define β_m≡ the particle speed in the plasma frame and μ_m≡ the cosine of the angle between the particle momentum and photon momentum in the plasma frame. Then

$$\begin{equation} \frac{dn_{m\gamma }}{dt} = \int \; d^3p \; \frac{dn_m}{d^3p} \; (1 - \mu _m \beta _m) \; \sigma c. \end{equation} \tag{ 26 }$$

It is convenient to rewrite this rate in terms of a "hot cross section"

$$\begin{equation} \sigma _h \equiv {1\over {n_m}} \int \; d^3p \; \frac{dn_m}{d^3p} \; (1 - \mu _m \beta _m) \; \sigma \end{equation} \tag{ 27 }$$

so that the interaction rate for a single photon is n_mσ_hc and the extinction coefficient is

$$\begin{equation} \alpha _{\nu } = n_m \sigma _h. \end{equation} \tag{ 28 }$$

So far we have assumed nothing about the interaction process.

5.1.2. Electron Scattering

For electron scattering the cross section is the Klein–Nishina total cross section expressed in terms of the photon energy in the electron rest frame ≡_e = γ_e(1 − μ_eβ_e), (γ_e ≡ (1 − β²_e)^−1/2 and we have substituted the subscript e for m):

$$\begin{equation} \sigma _{{\rm KN}}=\sigma _{T} \frac{3}{4 {\epsilon }_e^2} \left(2 + \frac{{\epsilon }_e^2 (1 + {\epsilon }_e)}{ (1 + 2{\epsilon }_e)^2} + \frac{{\epsilon }_e^2 - 2{\epsilon }_e - 2 }{2 {\epsilon }_e}\log (1+ 2{\epsilon }_e)\right). \end{equation} \tag{ 29 }$$

Here σ_T is the Thomson cross section. For ≪ 1,

$$\begin{equation} \sigma _{{\rm KN}} = \sigma _{T} (1 - 2 {\epsilon }+ O({\epsilon }^2)) \end{equation} \tag{ 30 }$$

which is numerically stable for small , unlike Equation (29).

Typically we assume a thermal electron distribution,

$$\begin{equation} \frac{dn_e}{d\gamma _e} = \frac{n_e}{\Theta _e}\frac{\gamma _e^2 \beta _e}{K_2\big(\Theta _e^{-1}\big)} \exp \left(-\frac{\gamma }{\Theta _e}\right), \end{equation} \tag{ 31 }$$

and evaluate (27) by direct integration to obtain σ_h(Θ_e, ). It is efficient to store the resulting cross sections in a two-dimensional lookup table at the beginning of the calculation. Our σ_h agrees with Wienke (1985).

5.2. Scattering Kernel

Once it is determined that a superphoton should be scattered at an event x^μ_s, the superphoton is passed to a scattering kernel which processes the scattering event according to the following procedure. Only unpolarized light is considered.

First, a plasma frame orthonormal tetrad is constructed by the same Gram–Schmidt orthogonalization procedure described in Section 2, and the old superphoton wave vector is transformed from the coordinate frame to the tetrad frame.

Second, a scattering electron is selected. We use the procedure described by Canfield et al. (1987), which selects the four-momentum p^μ_e of the scattering electron with an efficiency of 72% at Θ_e = 1 and nearly 100% for Θ_e ≪ 1 or Θ_e ≫ 1 (Canfield et al. 1987).

Third, we boost from the plasma frame tetrad to an electron frame tetrad q^μ_(a) and construct the scattered photon wave vector in this frame. The differential scattering cross section is sampled for the scattered photon energy '_e and the scattering angle θ. For low_energy photons (_e < _l; in grmonty  _l = 10⁻⁴) the scattering is approximately elastic, so we set '_e = _e and sample the Thomson differential cross section

$$\begin{equation} \frac{2\pi }{\sigma _T}\frac{d\sigma _T}{d\cos \theta } = \frac{3}{8} (1 + \cos ^2\theta) \end{equation} \tag{ 32 }$$

for the scattering angle θ using a rejection scheme. For _e>_l we sample the Klein–Nishina differential cross section

$$\begin{equation} \frac{2\pi }{\sigma _T} {{d\sigma _{{\rm KN}}}\over {d{\epsilon }_e^{\prime }}} = {1\over {{\epsilon }_e^2}}\left({{\epsilon }_e\over {{\epsilon }_e^{\prime }}} + {{\epsilon }_e^{\prime }\over {{\epsilon }_e}} - 1 + \cos ^2\theta\right), \end{equation} \tag{ 33 }$$

for '_e using a rejection scheme. Here cos θ = 1 + 1/_e − 1/'_e. This procedure is inefficient for _e ≫ 1, but in our target application such large photon energies are rare. Drawing a final angle ϕ from a uniform distribution on [0, 2π) completes the specification of the photon wave vector

$$\begin{equation} k^{(0)} = {\epsilon }_e^{\prime } \end{equation} \tag{ 34a }$$

$$\begin{equation} k^{(1)} = {\epsilon }_e^{\prime }\cos \theta \end{equation} \tag{ 34b }$$

$$\begin{equation} k^{(2)} = {\epsilon }_e^{\prime }\sin \theta \cos \phi \end{equation} \tag{ 34c }$$

$$\begin{equation} k^{(3)} = {\epsilon }_e^{\prime }\sin \theta \sin \phi \end{equation} \tag{ 34d }$$

in the electron-frame tetrad basis q^μ_(a); q₍₁₎ is aligned parallel to the spatial part of the incoming photon wave vector.

Finally, we boost back from the electron-frame tetrad to the plasma frame tetrad (some of these steps are combined in our code for computational efficiency), and use the plasma frame basis vectors to obtain the coordinate frame scattered photon wave vector k'_μ.

Spectra can be measured using a "detector" with area ΔA at distance R, frequency channels of logarithmic width Δln ν, and integration times ΔT. The flux density in frequency bin i is then just

$$\begin{equation} F_{\nu,i} = \frac{1}{\nu _i \Delta \ln \nu \; \Delta T\; \Delta A} \sum _j w_j (h \nu)_j, \end{equation} \tag{ 35 }$$

where the sum is taken over all superphotons j that land in the channel during the integration. In principle, a software detector can behave just like a physical detector, producing time-dependent spectra from time-dependent flow models. In practice, time-dependent models are not (yet) treated self-consistently.

To see why, consider a time-dependent model based on a general relativistic magnetohydrodynamics (GRMHD) model of a black hole accretion flow. Self-consistent treatment of the radiation field would require generating and tracking superphotons through the simulation data as it evolves, i.e., coupling the grmonty to the GRMHD code (the simulation data could be stored and post-processed, but this would require storing almost every time step and would be impractical and inefficient). The mean number of superphotons tracked simultaneously would depend on the desired signal to noise in the final spectrum, as well as the time and energy resolution. Our experience suggests that for nominal energy resolution, signal to noise, and time resolution of order the horizon light crossing time, ∼10⁸ superphotons would need to be maintained within the simulation for the duration of the evolution. The total number required is then N_tot ≫ 10⁸ and is currently inaccessible without significant computational resources. We plan to attempt this calculation later.

For now we construct spectra of time-dependent data using a stationary-data (or "fast-light") approximation: each time slice of data is treated as if it were stationary (time-independent). The slice emits superphotons for a time Δt. The photons then propagate through the slice data (as t varies along a geodesic the fluid variables are held fixed) and are detected at large distance. The time-steady flux is obtained by substituting Δt for ΔT in Equation (35).

6.1. Measuring νL_ν

In practice, we measured νL_ν,i rather than F_ν,i. Since νL_ν = 4πR²νF_ν, and R²/ΔA is the solid angle ΔΩ occupied by the detector,

$$\begin{equation} \nu L_{\nu,i} = {4\pi \over {\Delta \Omega \Delta t}} {1\over {\Delta \ln \nu }} \; \sum _j w_j h \nu _j. \end{equation} \tag{ 36 }$$

Typically our "detectors" capture all the superphotons in a large angular bin ΔΩ around the source. For example, in studies of axisymmetric black hole accretion flows the angular bins capture all photons with Boyer–Lindquist r>100GM/c², θ_n < θ < θ_n+1, independent of ϕ, where θ_n are the bin boundaries.

We can also estimate average values for any quantity Q associated with emission from a source (e.g., the absorption optical depth). We define the weight-averaged value of Q via

$$\begin{equation} {\langle } Q {\rangle } \equiv {{\sum Q w}\over {\sum w}}, \end{equation} \tag{ 37 }$$

where the sum is taken within an energy and angular bin.

6.2. Optimal Weights

The fractional variance in the Monte Carlo estimate for νL_ν is proportional to $\overline{w^2}/(\overline{N_s} \overline{w}^2)$, where overline means expectation value and N_s is the number of superphotons in the bin. Evidently the optimal weighting of superphotons is achieved when (1) the weights of superphotons are the same within bins and (2) the superphotons are evenly distributed across frequency bins.

If, as we created new superphotons, we knew νL_ν, then we could set

$$\begin{equation} w_\nu = \frac{1}{\overline{N_s} h} \; L_\nu \Delta \ln \nu \end{equation} \tag{ 38 }$$

and set $\overline{N_s} = N_{s,{\rm tot}}/N_b$, where N_b is the number of frequency bins.

Of course we do not know L_ν, but for the special case of an optically thin emitting plasma we can estimate it as

$$\begin{equation} L_{\nu } \approx \int d^3 x \sqrt{-g} \int d\Omega \; j_{\nu }. \end{equation} \tag{ 39 }$$

This estimate assumes that all photons escape to infinity, and it ignores Doppler shift, gravitational redshift, scattering, and angular structure in νL_ν. Nevertheless, it is useful because (1) it can be calculated before the Monte Carlo calculation begins, and (2) it is far better to use the information contained in this rough estimate of the spectrum than to proceed using, e.g., uniform weights.

We verify the accuracy of grmonty by comparing spectra produced on idealized problems against a reference spectrum (νL_ν)_ref computed analytically (when possible) or computed by an independent code. For all tests we use the following error norm, which effectively measures the maximum of the fractional error, compared to the reference solution, over frequency:

$$\begin{equation} {\langle \epsilon \rangle }= \frac{1}{\Delta \ln \nu } \int _{\nu _{{\rm min}}}^{\nu _{{\rm max}}} \frac{|(\nu L_\nu)_{{\rm grmonty}} - (\nu L_\nu)_{{\rm ref}}|}{(\nu L_\nu)_{{\rm ref}}} d\ln \nu, \end{equation} \tag{ 40 }$$

where Δln ν = ln(ν_max/ν_min) is the range of integration. The range of integration is the same as plotted in the spectra for each test.

7.1. Optically Thin Synchrotron Sphere

First we consider emission from a homogeneous, optically thin spherical cloud of unit volume threaded by a vertical magnetic field in flat space. The cloud parameters are Θ_e = 100, B = 1G, and n_e = 10¹⁵cm⁻³, which gives an optical depth at ν = 10⁹Hz of ∼10⁻² perpendicular to the magnetic field. The emissivity and absorptivity are constant along any line of sight, so

$$\begin{equation} I_{\nu } = \frac{j_{\nu }}{\alpha _{\nu }}(1 - e^{-\alpha _{\nu }L}) \approx j_{\nu } L + O\big(\tau _a^2\big), \end{equation} \tag{ 41 }$$

where L is the path length through the sphere. We numerically integrate this expression over detector solid angle to compute the spectral energy distribution, and compare with the spectrum grmonty produces in Figure 3. Evidently the result is unbiased. Figure 4 demonstrates that grmonty converges on the correct solution as ∝N^−1/2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

7.2. Optically Thick Synchrotron Sphere

The second test is identical except that the electron number density and therefore the optical depth are increased by a factor of 10⁵. The intensity along any line of sight is again

$$\begin{equation} I_{\nu } = \frac{j_{\nu }}{\alpha _{\nu }}(1 - e^{-\alpha _{\nu }L}). \end{equation} \tag{ 42 }$$

The emitting region becomes optically thick when (j_ν/B_ν)L ≳ 1. For the thermal synchrotron emissivity we use, this occurs below a critical frequency.

The spectrum is shown in Figure 5 and the convergence is shown in Figure 6. The figures make two key points: convergence is slow for small numbers of superphotons; and the overall magnitude of the error is larger than in the optically thin case shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The slow initial convergence is due to the large optical depth at some frequencies. When the optical depth is large no superphotons of appreciable weight are recorded until some superphotons have been created in the fraction ∼1/τ of the volume that lies within the photosphere. Our problem has τ ∼ 10⁵ at ν ∼ 10⁸Hz. Since 〈〉 effectively measures the maximum of the error over frequency, it is not surprising that grmonty requires ∼10⁶ superphotons before it begins to converge as N^−1/2.

7.3. Comptonization of Soft Photons in a Spherical Cloud of Plasma

This test is based on a problem posed in Section 6 of Pozdynakov et al. (1983): the spectrum of a spherical, homogeneous, unmagnetized cloud of radius R that contains thermal electrons at density n_e and temperature T_e that scatter light from a central, thermal source of temperature T_s. Absorption and emission in the cloud are neglected. The dimensionless parameters of the problem are Θ_e = kT_e/(m_ec²), τ = Rσ_Tn_e, and Θ_s = kT_s/(m_ec²).

For this test our reference spectrum is computed with an implementation of Pozdnyakov et al.'s Monte Carlo scheme kindly provided by S. Davis. This code, sphere, has been modified in two ways: we have replaced the approximate hot cross sections defined in Pozdynakov et al. (1983) with our more exact, numerically integrated values, and we use the exact Klein–Nishina cross section Equation (29) when choosing the electron with which a photon should scatter. Without these changes to sphere differences between the spectra are ≲1%, consistent with the error Pozdynakov et al. (1983) quote for their approximations. These small differences are enough, however, to prevent grmonty from converging as expected for large numbers of superphotons.

Figures 7–9 show the radiation spectra (upper panels) produced by both codes and a fractional difference (bottom panels) between them for Θ_e = 4, Θ_s = 10⁻⁸ and for various values of the optical depth τ = 10⁻⁴, 0.1, and 3. Figure 10 demonstrates that grmonty converges to the reference solution as ∝N^−1/2 for each optical depth.

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

7.4. Synchrotron Self-absorbed Spectra in Black Hole Spacetimes

We consider two idealized problems: (1) smooth, spherically symmetric infall onto a Schwarzschild black hole; (2) a snapshot of a turbulent accretion flow around a Kerr black hole with a/M = 0.9375 produced by general relativistic MHD simulation with the HARM code (Gammie et al. 2003). In this test, our reference solution is computed by the ray tracing code ibothros (Noble et al. 2007). ibothros solves the invariant form of the radiative transfer equation along geodesics that terminate in a fictitious camera at large distance. Spectra are constructed by imaging the source at successive frequencies and performing an angular integral over the images to estimate νL_ν. In these tests, scattering is turned off in grmonty.

There are algorithmic differences between grmonty and ibothros that lead to differences in their spectra.

First, grmonty measures the flux in energy and angular bins, whereas ibothros measures the flux at a particular inclination and energy. This is not an important effect unless there is sharp angular or energy structure in the spectrum.

The next difference is more subtle and is related to the treatment of gridded model data used to construct these tests. In grmonty,  quantities such as the density, temperature, etc., are viewed as the average of these variables over a grid zone. In ibothros,  the grid variables are viewed as zone-centered samples and a continuous distribution is created by multi-linear interpolation between zone centers. The difference is illustrated in one dimension in Figure 11. ibothros and grmonty therefore differ in zone-averaged emissivity by O(Δx/L)², where Δx is the zone size and L is the characteristic scale of the emitting structure.

Zoom In Zoom Out Reset image size

Differences in the grmonty and ibothros spectra of structures with Δx ≪ L are therefore small. High-frequency synchrotron emission, however, is exponentially dominated by emission from a few zones with highest ν_s (see Equation (4b); these are the zones with highest temperature or strongest magnetic field). Then Δx/L ∼ 1 for high-frequency emission, and the grmonty and ibothros spectra differ by of order unity. A similar effect occurs for low-frequency synchrotron emission where the optical depth is large. The spectrum is sensitive to the run of physical variables through the photosphere, which has size comparable to or smaller than a grid zone. The subsequent test results will omit these parts of the spectra. These inconsistencies between grmonty and ibothros could be eliminated by using identical, continuous models for subgrid reconstruction. But this would require significant investment in recoding that, in our view, is not worthwhile: to the extent that the spectrum depends on the flow structure at and below the grid scale, it is not reliable.

In spite of these differences, and other differences between the codes related to differing accuracy parameters, grmonty should converge on the ibothros result until the effects of data interpolation become the dominant sources of error.

Figures 12 and 13 show the spectrum and convergence relative to ibothros, respectively, for the spherical accretion problem. This problem is an attractive test for at least two reasons. First, the emission is isotropic so the effects of angular binning in grmonty are eliminated. Second, the flow is smooth, so the differences due to data interpolation will be small. Evidently the spherical accretion model converges at the expected, N^−1/2, rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 14 and 15 show the spectrum and convergence relative to ibothros, respectively, for the turbulent accretion problem. This comparison is more challenging because the high-energy emission originates in compact "hot spots." Evidently the two models agree at the 10% level everywhere, with the largest differences at high and low frequencies, where the subgrid reconstruction comes into play. Excluding these high- and low-frequency regions, the agreement between the codes is at the few percent level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now apply grmonty to the same HARM simulation data used in the grmonty-ibothros comparison above, this time with scattering enabled. Figure 16 shows the resulting spectrum with the ibothros result shown for comparison. Evidently, scattering has little effect on the sub-mm spectrum since the scattering optical depth is small (∼10⁻⁴) but the model now predicts a significant X-ray flux. A more detailed analysis of Comptonized spectra from GRMHD simulations in the context of Sgr A* will be given in a separate paper (M. Mościbrodzka et al. 2009, in preparation).

Zoom In Zoom Out Reset image size

Figure 17 shows the results of self-convergence tests for the full code. Convergence is initially slow but quickly approaches a rate proportional to N^−1/2 as spectral bins become sufficiently sampled. The spectra shown in Figure 16, which were taken from a single run with 10⁹ superphotons, were used for references.

Zoom In Zoom Out Reset image size

We have described and tested a code that solves the radiative transfer problem for optically thin ionized plasmas in general spacetimes. The code treats the full angular dependence of emission and absorption, treats single Compton scattering exactly (double Compton and induced Compton scattering are neglected), and can be easily adapted to simulate emission from both analytic and numerical models. While we have specialized to synchrotron emission in this work, grmonty is constructed so that it is straightforward to include other relevant emission mechanisms such as bremsstrahlung with only minimal modification.

As a demonstration of a practical use for grmonty, we have computed the first spectra, including synchrotron emission and Compton scattering, from GRMHD models of a turbulent accretion disk. Other potential applications of our code are to neutron star accretion, emission from relativistic blast waves, and any problem where relativistic bulk motion makes radiation transfer treatments that expand the flow in orders of v/c problematic.

This work was supported by the National Science Foundation under grants AST 00-93091, PHY 02-05155, and AST 07-09246, and by a Richard and Margaret Romano Professorial scholarship, a Sony faculty fellowship, and a University Scholar appointment to C.F.G. Portions of this work were performed while C.F.G. was a Member at the Institute for Advanced Study in academic year 2006–2007. The authors are especially grateful to Shane Davis and Stu Shapiro for their insights, and to Peter Goldreich, Fred Lamb, and John Hawley for discussions. We thank the anonymous referee for suggestions which helped improve this manuscript in clarity and completeness.