Blacklight: A General-relativistic Ray-tracing and Analysis Tool

There are two common types of camera used in general-relativistic ray tracing, which we will term "plane-parallel" and "pinhole." The former takes pixels to be located at distinct points in a plane (suitably defined), each seeing light with the same parallel (again, suitably defined) momentum. The latter takes all pixels to be collocated in space but sensitive to different momentum directions. Plane-parallel cameras are appropriate for making images as seen at infinity without having to place the camera at large distances from the source, while pinhole cameras model what a small detector would see at its location. Codes in the literature use one or the other, and Blacklight supports both.

In either case, Blacklight requires the user to specify the position xⁱ of the camera center, the normal-frame velocity u^p of the camera center, and the unnormalized momentum ${k}_{i}={\delta }_{i}^{p}{k}_{p}$ of the light received at the camera center. Velocities are specified as u^p rather than uⁱ , as the latter fails to uniquely determine a future-directed 4-velocity within the ergosphere.

With the components of x, u, and k known in any coordinate system, we can define unit vectors for the line-of-sight, vertical, and horizontal directions. To do this, we work in the camera frame based on Cartesian Kerr–Schild coordinates, where the metric has components

$$\begin{eqnarray}&&{g}_{t^{\prime} a^{\prime} }=-{\delta }_{{ta}},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{g}_{a^{\prime} b^{\prime} }={g}_{{ab}}-\displaystyle \frac{{u}_{a}}{{u}_{t}}{g}_{{tb}}-\displaystyle \frac{{u}_{b}}{{u}_{t}}{g}_{{ta}}+\displaystyle \frac{{u}_{a}{u}_{b}}{{u}_{t}{u}_{t}}{g}_{{tt}},\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&{g}^{t^{\prime} a^{\prime} }=-{\delta }^{{ta}},\end{eqnarray} \tag{ 7c }$$

$$\begin{eqnarray}&&{g}^{a^{\prime} b^{\prime} }={g}^{{ab}}+{u}^{a}{u}^{b}.\end{eqnarray} \tag{ 7d }$$

The line-of-sight vector K is a unit vector parallel to k. Explicitly, we have

$$\begin{eqnarray}&&{K}_{a^{\prime} }\propto {k}_{a}-\displaystyle \frac{{u}_{a}}{{u}_{t}}{k}_{t},\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{K}^{t^{\prime} }\propto -{u}^{\alpha }{k}_{\alpha },\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&{K}^{a^{\prime} }\propto {g}^{a^{\prime} b^{\prime} }{K}_{b^{\prime} },\end{eqnarray} \tag{ 8c }$$

where the proportionality constant is set by requiring ${K}_{a^{\prime} }{K}^{a^{\prime} }=1$. The transformation rule

$$\begin{eqnarray}&&{K}^{t}={u}^{t}{K}^{t^{\prime} }-\displaystyle \frac{{u}_{a}}{{u}_{t}}{K}^{a^{\prime} },\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&{K}^{a}={K}^{a^{\prime} }+{u}^{a}{K}^{t^{\prime} }\end{eqnarray} \tag{ 9b }$$

returns these components to the coordinate frame.

A vertical direction v begins with a fiducial "up" vector ${U}^{a^{\prime} }={\delta }_{z}^{a}$. In the special case that the camera is centered on the polar (z) axis, we instead take ${U}^{a^{\prime} }={\delta }_{y}^{a}$. Then v is defined from U orthogonal to K via Gram–Schmidt,

$$\begin{eqnarray}&&{v}^{t^{\prime} }\propto 0,\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{v}^{a^{\prime} }\propto {U}^{a^{\prime} }-{U}^{b^{\prime} }{K}_{b^{\prime} }{K}^{a^{\prime} },\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&{v}_{a^{\prime} }\propto {g}_{a^{\prime} b^{\prime} }{v}^{b^{\prime} },\end{eqnarray} \tag{ 10c }$$

where the normalization is determined by ${v}_{a^{\prime} }{v}^{a^{\prime} }=1$. The orthogonal horizontal direction is defined as

$$\begin{eqnarray}&&{h}^{t^{\prime} }=0,\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&{h}^{a^{\prime} }=\displaystyle \frac{1}{\sqrt{D}}[a\ b\ c]{v}_{b^{\prime} }{K}_{c^{\prime} },\end{eqnarray} \tag{ 11b }$$

where D is the determinant of ${g}_{a^{\prime} b^{\prime} }$ and [a b c] is the Levi–Civita symbol. Blacklight allows users to rotate the camera about its axis by an angle ψ, so that the horizontal and vertical directions become

$$\begin{eqnarray}&&H=h\cos \psi -v\sin \psi ,\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&V=v\cos \psi +h\sin \psi .\end{eqnarray} \tag{ 12b }$$

Users specify a camera width w and linear resolution N. Consider a pixel with horizontal and vertical indices i, j ∈{0,...,N − 1}. Define the coefficients

$$\begin{eqnarray}&&a,b=\left(\{i,j\}-\displaystyle \frac{N}{2}+\displaystyle \frac{1}{2}\right)\displaystyle \frac{w}{N}.\end{eqnarray} \tag{ 13 }$$

For plane-parallel cameras, the pixel position and momentum are defined to be

$$\begin{eqnarray}&&{x}_{\mathrm{pix}}^{\alpha }={x}^{\alpha }+{\rm{\Delta }}{x}^{\alpha },\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{k}_{\mathrm{pix}}^{a}={K}^{a},\end{eqnarray} \tag{ 14b }$$

where the displacement in the camera frame is

$$\begin{eqnarray}&&{\rm{\Delta }}{x}^{\alpha ^{\prime} }={{aH}}^{\alpha ^{\prime} }+{{bV}}^{\alpha ^{\prime} }.\end{eqnarray} \tag{ 15 }$$

For pinhole cameras, they are

$$\begin{eqnarray}&&{x}_{\mathrm{pix}}^{\alpha }={x}^{\alpha },\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&{k}_{\mathrm{pix}}^{a}={k}_{\mathrm{pix}}^{a^{\prime} }+{u}^{a}{k}_{\mathrm{pix}}^{t^{\prime} },\end{eqnarray} \tag{ 16b }$$

where we have

$$\begin{eqnarray}&&{k}_{\mathrm{pix}}^{t^{\prime} }={K}^{t},\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{k}_{\mathrm{pix}}^{a^{\prime} }=\displaystyle \frac{{{rK}}^{a}-{{aH}}^{a}-{{bV}}^{a}}{{({r}^{2}+{a}^{2}+{b}^{2})}^{1/2}}.\end{eqnarray} \tag{ 17b }$$

For both cameras, the remaining component ${k}_{\mathrm{pix}}^{t}$ is found such that ${g}_{\alpha \beta }{k}_{\mathrm{pix}}^{\alpha }{k}_{\mathrm{pix}}^{\beta }=0$. There is a unique positive root on and outside the ergosphere. Inside the ergosphere, the lesser of the two positive roots is chosen. Note the normalization of k_pix is done in an arbitrary but well-defined way here without reference to the frequency of light being observed, fixing the scale of the affine parameter used in the geodesic equation.

Eventually, the physical frequency ν (e.g., 230 GHz) must come into play. At the camera, the ray's physical momentum will be nk_pix, where the normalization n is set by one of two user options. If the camera is to record frequency ν at its location (most appropriate for pinhole cameras), then

$$\begin{eqnarray}&&-{{nk}}_{\alpha }^{\mathrm{pix}}{u}^{\alpha }=\nu .\end{eqnarray} \tag{ 18 }$$

Alternatively, if the camera is to record light that would be observed at frequency ν by a static observer at infinity (most appropriate for plane-parallel cameras), then

$$\begin{eqnarray}&&-{{nk}}_{t}^{\mathrm{pix}}=\nu .\end{eqnarray} \tag{ 19 }$$

Given initial contravariant position components x^α and covariant null momentum components k_α , Blacklight integrates the geodesic equation in the Hamiltonian form

$$\begin{eqnarray}&&\displaystyle \frac{{d}{x}^{\alpha }}{{d}\lambda }={g}^{\alpha \beta }{k}_{\beta },\end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{d}{k}_{t}}{{d}\lambda }=0,\end{eqnarray} \tag{ 20b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{d}{k}_{a}}{{d}\lambda }=-\displaystyle \frac{1}{2}{\partial }_{a}{g}^{\alpha \beta }{k}_{\alpha }{k}_{\beta }.\end{eqnarray} \tag{ 20c }$$

(see James et al. 2015; Pu et al. 2016; Kawashima et al. 2021; Velásquez-Cadavid et al. 2022). It also sometimes tracks proper distance along the ray via

$$\begin{eqnarray}&&\displaystyle \frac{{d}s}{{d}\lambda }={\left({g}_{{pq}}{k}^{p}{k}^{q}\right)}^{1/2}.\end{eqnarray} \tag{ 21 }$$

Integration proceeds in the direction of negative λ (backward in time) from the camera. It terminates when the radial coordinate r crosses the threshold from less than to greater than that of the camera center. Even in horizon-penetrating coordinates, backward-propagating rays cannot reach the horizon in finite affine parameter. Geodesics are therefore also terminated if they cross inside an inner r threshold set by the user.

Blacklight offers three numerical schemes for integrating Equation (20a). One is the standard second-order Runge–Kutta (RK2) Heun method, which has the Butcher tableau given in Table 1, and another is the standard fourth-order Runge–Kutta (RK4) method with the Butcher tableau of Table 2. In both cases, the step size in affine parameter is taken to be

$$\begin{eqnarray}&&{\rm{\Delta }}\lambda ={f}_{\mathrm{step}}(r-{r}_{\mathrm{hor}}),\end{eqnarray} \tag{ 22 }$$

where f_step is a user-specified constant and r_hor is the outer horizon radius.

Table 1. Butcher Tableau for the RK2 Geodesic Integrator

Download table as:  ASCIITypeset image

Table 2. Butcher Tableau for the RK4 Geodesic Integrator

Download table as:  ASCIITypeset image

In terms of robustness and ease of use, however, the recommended geodesic integrator is the "RK5(4)7M" method of Dormand & Prince (1980; DP), which is an adaptive fifth-order Runge–Kutta method that calculates RK4 values from the same set of function evaluations in order to estimate errors for each step. In comparing common integrators in the context of the computation of null geodesics around black holes, Velásquez-Cadavid et al. (2022) find DP to be far more accurate. The step size is adjusted dynamically based on this error estimate, in order to take steps as large as possible given a user-prescribed tolerance (see Press et al. 1997). Here the error is defined in terms of the initial, fourth-order, and fifth-order values of the eight dependent variables {y_m } = {x^α } ∪ {k_α } in (20a) via

$$\begin{eqnarray}&&\epsilon =\mathop{\max }\limits_{m}\left\{\displaystyle \frac{\left|{y}_{m}^{(5)}-{y}_{m}^{(4)}\right|}{{f}_{\mathrm{abs}}+{f}_{\mathrm{rel}}\max \left(\left|{y}_{m}^{(\mathrm{init})}\right|,\left|{y}_{m}^{(5)}\right|\right)}\right\},\end{eqnarray} \tag{ 23 }$$

where f_abs and f_rel are provided by the user. Values of > 1 indicate a step must be reattempted with a smaller Δλ.

As described in Shampine (1986), the function evaluations can be arranged in a different way to provide a fourth-order accurate estimate of the midpoint of a step. Blacklight takes this midpoint to be the position and momentum of the step for the purposes of radiative transfer integration. In cases where the geodesic steps cover more proper distance than a user-specified limit (such that the radiative transfer equation would be poorly sampled even where the geodesic is calculated accurately), Blacklight subdivides the step. The aforementioned midpoint, together with function values and derivatives at both endpoints of the step, can be used to uniquely define a quartic expression, providing a fourth-order accurate interpolation anywhere within the step (Shampine 1986) and enabling this subdivision. The exact choice of coefficients used in DP implementations varies; the ones used here are given in Table 3.

Table 3. Butcher Tableau for the DP Geodesic Integrator

Download table as:  ASCIITypeset image

In all cases, integration is performed in Cartesian Kerr–Schild coordinates. As these coordinates are stationary, use of Equation (20a) ensures the covariant energy is conserved along rays exactly in the code, as it is physically. Because Blacklight does not employ spherical coordinates here, exact conservation of ϕ covariant momentum (z angular momentum) is not guaranteed. However, this choice avoids numerical issues arising from geodesics approaching coordinate poles.

Blacklight can create unpolarized images, ignoring polarization effects along each ray. In terms of the Lorentz invariants

$$\begin{eqnarray}&&I={I}_{\nu }{\nu }^{-3},\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&{j}_{I}={j}_{\nu }{\nu }^{-2},\end{eqnarray} \tag{ 24b }$$

$$\begin{eqnarray}&&{\alpha }_{I}={\alpha }_{\nu }\nu ,\end{eqnarray} \tag{ 24c }$$

the general-relativistic equation for unpolarized radiative transfer is

$$\begin{eqnarray}&&\displaystyle \frac{{d}I}{{d}\lambda }={j}_{I}-{\alpha }_{I}I,\end{eqnarray} \tag{ 25 }$$

where the affine parameter λ must be normalized such that Equation (20a) holds. We take j_I and α_I to be constant over each ray segment (typically no larger than nearby cells in the underlying simulation). The evolution of I from the beginning of a segment (where we denote it I⁻) to the end (where we denote it I⁺), across an optical depth Δτ = α_I Δλ, employs the well-known exact solution

$$\begin{eqnarray}{I}^{+}=\left\{\begin{array}{ll}\exp (-{\rm{\Delta }}\tau )\left({I}^{-}+\mathrm{expm}1({\rm{\Delta }}\tau )\displaystyle \frac{{{\rm{j}}}_{I}}{{\alpha }_{I}}\right), & {\alpha }_{I}\gt 0\ \mathrm{and}\ {\rm{\Delta }}\tau \lt {\rm{\Delta }}{\tau }_{\max };\\ \displaystyle \frac{{j}_{I}}{{\alpha }_{I}}, & {\alpha }_{I}\gt 0\ \mathrm{and}\ {\rm{\Delta }}\tau \geqslant {\rm{\Delta }}{\tau }_{\max };\\ {I}^{-}+{j}_{I}{\rm{\Delta }}\lambda , & {\alpha }_{I}=0;\end{array}\right.\end{eqnarray} \tag{ 26 }$$

where we use the standard function $\mathrm{expm}1(\cdot )=\exp (\cdot )-1$ as written and fix a finite ${\rm{\Delta }}{\tau }_{\max }=100$ to avoid floating-point inaccuracies.

Users can choose to integrate the equations of polarized radiative transfer, evolving four degrees of freedom and producing images in Stokes parameters I_ν , Q_ν , U_ν , and V_ν at the camera. In addition to j_I and α_I , this requires knowing the Lorentz-invariant linearly polarized emissivities j_Q and j_U , circularly polarized emissivity j_V , linearly polarized absorptivities α_Q and α_U , circularly polarized absorptivity α_V , Faraday conversion rotativities ρ_Q and ρ_U , and Faraday rotation rotativity ρ_V .

There are two primary approaches used by polarized general-relativistic ray-tracing codes. Some (e.g., grtrans) make use of the convenient form of the Walker–Penrose constant (Walker & Penrose 1970) along null geodesics in Boyer–Lindquist coordinates to parallel-transport fiducial polarization directions backward from the camera, thus allowing each segment's local definition of the Stokes parameters to match those at the camera. Here we adopt the alternative approach of evolving the polarized state of the radiation forward in a way commensurate with both plasma and coordinate effects, as done in, e.g., ipole. This allows our use of Cartesian (or even more general) coordinates.

The method we use is well described by Mościbrodzka & Gammie (2018); we only briefly summarize it here. The polarized transfer equation can be written as both

$$\begin{eqnarray}&&\displaystyle \frac{{d}{S}_{A}}{{d}\lambda }=\mathop{\underbrace{{J}_{A}-{M}_{{AB}}{S}_{B}}}\limits_{\mathrm{plasma}\ {\rm{terms}}}\ +\ (\mathrm{coordinate}\ {\rm{terms}})\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{d}{N}^{\alpha \beta }}{{d}\lambda }=\mathop{\underbrace{-{k}^{\gamma }\left({{\rm{\Gamma }}}_{\gamma \delta }^{\alpha }{N}^{\delta \beta }-{{\rm{\Gamma }}}_{\gamma \delta }^{\beta }{N}^{\alpha \delta }\right)}}\limits_{\mathrm{coordinate}\ {\rm{terms}}}\ +\ (\mathrm{plasma}\ {\rm{terms}}).\end{eqnarray} \tag{ 28 }$$

Here A and B run over the Stokes indices (with repetition implying summation), S is the vector of Stokes parameters, J is the vector of emissivities, M is the matrix of absorptivities and rotativities, N is the complex Hermitian coherency tensor, and Γ is the collection of connection coefficients. The dependent variables are encoded in both S and N, either of which can be used to define the other in a given frame.

The general procedure is to use Strang splitting for each ray segment of affine parameter length Δλ, evolving Equation (28) without the unwieldy plasma terms for Δλ/2, then evolving Equation (27) without the equally unwieldy coordinate terms for Δλ, and finally evolving Equation (28) without plasma coupling for another Δλ/2. All coefficients are held constant over each substep. As the last substep in one segment and the first substep in the next segment evolve the same equation, we merge these two via a second-order predictor-corrector method with the Butcher tableau given in Table 4.

Table 4. Butcher Tableau for the Predictor-corrector Evolution of Equation (28)

Download table as:  ASCIITypeset image

Just as Equation (25) requires special care and exact solutions to avoid numerically unstable integration, so too does Equation (27). The full solution is described in a number of sources, such as Jones & O'Dell (1977), Landi Degl'Innocenti & Landi Degl'Innocenti (1985), and Mościbrodzka & Gammie (2018). ⁷ The latter provides useful simplifying cases for when various coefficients vanish, which Blacklight employs when applicable.

Given a simulation grid, possibly including mesh refinement, and a set of rays with sample points, Blacklight interpolates values from the former onto the latter. This interpolation is done in the primitive variables employed by Athena++: total rest-mass density ρ; total pressure p or electron entropy κ_e, depending on electron model (see Section 2.6); fluid velocity u^p , or the Cartesian Kerr–Schild normal-frame values if appropriate; and magnetic field Bⁱ or B^a , as appropriate.

Interpolation can be done in nearest-neighbor fashion, where the values are copied from whatever unique cell bounds the sample point. Alternatively, Blacklight can perform trilinear interpolation in simulation coordinates. Athena++ domain-decomposes the simulation into logically Cartesian blocks of cells, and interpolation within a block costs negligibly more than nearest-neighbor lookup. Accurate interpolation across block boundaries, however, can dominate the cost of producing an image, given that adjacent blocks may be at different refinement levels and ghost zones are generally not included in data dumps. Interblock interpolation rarely makes a noticeable difference in the results, compared to intrablock interpolation; users can choose which method to employ.

The plasma state at each sample point is fully defined after interpolation. Useful quantities that can be derived from the interpolated ones, possibly with user-defined assumptions about the plasma, include electron number density n_e; electron dimensionless temperature Θ_e = k_B T_e/m_e c², with T_e the electron temperature; fluid-frame magnetic field magnitude $B={({b}_{\alpha }{b}^{\alpha })}^{1/2}$, with b^α = u_β (*F)^{β α} the components of the magnetic 4-vector and F the electromagnetic field tensor; plasma σ = B²/ρ; and plasma β⁻¹ = B²/2p.

It is standard practice to rotate one's coordinate system at each point along a ray in order to make j_U , α_U , and ρ_U vanish. For the remaining eight coefficients, Blacklight uses the fitting formulas provided by Marszewski et al. (2021), which conveniently unifies and corrects the synchrotron literature. These formulas cover the cases of thermal, power-law, and kappa distributions for the electrons, always assumed to be isotropic. Blacklight can calculate the coefficients resulting from any admixture of these three populations. Thermal electrons are distributed in Lorentz factor γ according to the Maxwell–Jüttner distribution based on Θ_e,

$$\begin{eqnarray}&&\displaystyle \frac{{d}{n}_{{\rm{e}}}}{{d}\gamma }=\displaystyle \frac{{n}_{{\rm{e}}}\gamma \sqrt{{\gamma }^{2}-1}}{{{\rm{\Theta }}}_{{\rm{e}}}{K}_{2}({{\rm{\Theta }}}_{{\rm{e}}}^{-1})}\exp \left(-\displaystyle \frac{\gamma }{{{\rm{\Theta }}}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 29 }$$

where K₂ is the cylindrical Bessel function. Power-law electrons are defined by an index p and Lorentz-factor cutoffs ${\gamma }_{\min }$ and ${\gamma }_{\max }$,

$$\begin{eqnarray}\displaystyle \frac{{d}{n}_{{\rm{e}}}}{{d}\gamma }\propto \left\{\begin{array}{ll}{\gamma }^{-p}, & \gamma \geqslant {\gamma }_{\min }\ \mathrm{and}\ \gamma \leqslant {\gamma }_{\max };\\ 0, & \gamma \lt {\gamma }_{\min }\ \mathrm{or}\ \gamma \gt {\gamma }_{\max }.\end{array}\right.\end{eqnarray} \tag{ 30 }$$

The kappa distribution is defined by index κ and width w,

$$\begin{eqnarray}&&\displaystyle \frac{{d}{n}_{{\rm{e}}}}{{d}\gamma }\propto \gamma \sqrt{{\gamma }^{2}-1}{\left(1+\displaystyle \frac{\gamma -1}{\kappa w}\right)}^{-\kappa -1}.\end{eqnarray} \tag{ 31 }$$

In all cases, the electron number density is given by

$$\begin{eqnarray}&&{n}_{{\rm{e}}}=\displaystyle \frac{\rho }{\mu {m}_{{\rm{p}}}}{\left(1+\displaystyle \frac{1}{{n}_{{\rm{e}}}/{n}_{{\rm{i}}}}\right)}^{-1},\end{eqnarray} \tag{ 32 }$$

where the user specifies the molecular weight μ and electron-to-ion number-density ratio n_e/n_i.

Unlike the power-law and kappa-distribution parameters, T_e varies across the simulation domain. As the simulations for which Blacklight is intended often treat the plasma as a single fluid with a single temperature T, we must invoke a prescription for determining T_e as a function of the local plasma state available in postprocessing. The code supports the common choice of declaring a plasma-beta-dependent ratio between the ion temperature and electron temperature,

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{i}}}}{{T}_{{\rm{e}}}}=\displaystyle \frac{{R}_{\mathrm{high}}+{\beta }^{-2}{R}_{\mathrm{low}}}{1+{\beta }^{-2}}\end{eqnarray} \tag{ 33 }$$

with global constants R_high and R_low, as introduced by Mościbrodzka et al. (2016). In this case, the electron temperature is

$$\begin{eqnarray}&&{T}_{{\rm{e}}}=\displaystyle \frac{\mu {m}_{{\rm{p}}}p}{{k}_{{\rm{B}}}\rho }\cdot \displaystyle \frac{{n}_{{\rm{e}}}/{n}_{{\rm{i}}}+1}{{n}_{{\rm{e}}}/{n}_{{\rm{i}}}+{T}_{{\rm{i}}}/{T}_{{\rm{e}}}}.\end{eqnarray} \tag{ 34 }$$

Alternatively, Blacklight can use the electron entropy available in two-temperature simulations as prescribed by Ressler et al. (2015). In this case, the code must provide the values of an entropy-like variable ${\kappa }_{{\rm{e}}}\propto \exp ({s}_{{\rm{e}}}/{k}_{{\rm{B}}})$, where s_e is the electron entropy per electron. Then the electron temperature is given by

$$\begin{eqnarray}&&{T}_{{\rm{e}}}=\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{5{k}_{{\rm{B}}}}\left({\left(1+25{({m}_{{\rm{e}}}{n}_{{\rm{e}}}{\kappa }_{{\rm{e}}})}^{2/3}\right)}^{1/2}-1\right)\end{eqnarray} \tag{ 35 }$$

(see Sądowski et al. 2017).

For null geodesics confined to the midplane in Kerr spacetime and outside the photon ring(s), solutions to the geodesic equation involve relatively tractable integrals. In Gold et al. (2020), the codes' values of the azimuthal deflection angles Δϕ of such geodesics are compared to the values given by formulas in Iyer & Hansen (2009) for the asymptotic deflection from past to future infinity. Here we repeat this test with the following modification. Given that the geodesic origin (source) and termination point (camera) are at finite distance (r = 1000 GM/c² here), we compare the result of ray tracing to the finite-distance formulas from Iyer & Hansen,

$$\begin{eqnarray}&&{\rm{\Delta }}\phi ={\int }_{{u}_{\mathrm{source}}}^{{u}_{0}}\left|\displaystyle \frac{{d}\phi }{{d}u}\right|\,{d}u-{\int }_{{u}_{0}}^{{u}_{\mathrm{camera}}}\left|\displaystyle \frac{{d}\phi }{{d}u}\right|\,{d}u,\end{eqnarray} \tag{ 36a }$$

$$\begin{eqnarray}\begin{array}{rcl}\left|\displaystyle \frac{{d}\phi }{{d}u}\right| & = & \displaystyle \frac{1-2M(1-a/b)u}{1-2{Mu}+{a}^{2}{u}^{2}}\\ & & \times \,{\left(2M{\left(1-\displaystyle \frac{a}{b}\right)}^{2}{u}^{3}-\left(1-\displaystyle \frac{{a}^{2}}{{b}^{2}}\right){u}^{2}+\displaystyle \frac{1}{{b}^{2}}\right)}^{-1/2},\end{array}\end{eqnarray} \tag{ 36b }$$

$$\begin{eqnarray}\begin{array}{rcl}{r}_{0} & = & \displaystyle \frac{1}{{u}_{0}}=\displaystyle \frac{2\left|b\right|}{\sqrt{3}}{\left(1-\displaystyle \frac{{a}^{2}}{{b}^{2}}\right)}^{1/2}\\ & & \times \,\cos \left(\displaystyle \frac{1}{3}{\cos }^{-1}\left(-\displaystyle \frac{3\sqrt{3}M{(1-a/b)}^{2}}{\left|b\right|{(1-{a}^{2}/{b}^{2})}^{3/2}}\right)\right),\end{array}\end{eqnarray} \tag{ 36c }$$

where u = 1/r, r₀ is the radial coordinate of closest approach, and b = − k_ϕ /k_t is the impact parameter.

We place a camera centered at r = 1000 GM/c², θ = π/2, and ϕ = 0 around a black hole with a/M = 0.9. There are 51 pixels across the middle of the image, spanning a width w = 36 GM/c². Of these pixels, 14 are ignored in this analysis; they lie inside the photon ring, given they have impact parameters satisfying b₋ ≤ b ≤ b₊ for the critical values

$$\begin{eqnarray}&&{b}_{\pm }=-a\pm 6M\cos \left(\displaystyle \frac{1}{3}{\cos }^{-1}\left(\mp \displaystyle \frac{a}{M}\right)\right).\end{eqnarray} \tag{ 37 }$$

Figure 1 shows the deflection angles calculated by Blacklight compared to Equation (36). The left and right sides of the figure correspond to the left and right sides of the image. For plotting purposes, we subtract the asymptotic, flat-spacetime "deflection" of ${\rm{\Delta }}{\phi }_{\mathrm{flat}}=\mathrm{sgn}({\rm{b}})\pi$. Here we employ the DP integrator with tolerances f_abs, f_rel = 10⁻⁸, and the agreement is very good, with absolute errors in the range 10⁻⁶–10⁻⁵. A similar test in Gold et al. (2020, see their Figure 1) shows most codes' errors being 10⁻⁵ or greater, though in that case, the finite-distance code values are compared to the asymptotic formula.

Zoom In Zoom Out Reset image size

We can compare the errors obtained with other numerical tolerances and integration schemes. The top panel of Figure 2 shows the errors for the DP integrator using tolerances ranging over 10⁻¹²–10⁻⁵, plotted against the number of steps taken to trace each ray. The adaptive nature of this integrator means stricter tolerances can be satisfied by only using slight more points near closest approach. For this reason, we recommend DP over the RK integrators for general-purpose ray tracing, where reasonable tolerances set by the user almost always yield reasonably accurate results without the danger of requiring intractably many steps.

Zoom In Zoom Out Reset image size

When testing the RK4 and RK2 integrators, we vary f_step over 0.002–0.1. The solution improves with decreasing step size, as expected, sometimes even outperforming typical DP results. However, this accuracy can come at a large cost in terms of the number of steps taken. Conversely, the RK integrators allow very quick and somewhat inaccurate integrations using relatively few steps. Relative to DP, they permit users to make larger trade-offs between accuracy and performance. It is possible that these non-adaptive variable-step-size methods could be improved for this type of problem with a better heuristic than Equation (22).

Gold et al. (2020) compare the images produced by unpolarized transfer using five sets of parameterized formulas for velocity, emissivity, and absorptivity, shown in their Figure 2. We repeat these tests with the same parameters and camera settings as in that work, choosing the DP integrator, a plane-parallel camera, and frequencies of 230 GHz at the location of the camera. The resulting images, shown in Figure 3, agree with those obtained by other codes.

Zoom In Zoom Out Reset image size

We report the total fluxes for these images in Table 2. Comparing to Table 2 of Gold et al. (2020), Blacklight lies comfortably in the middles of the ranges spanned by the seven codes in that work. Our table shows how the minimum and maximum values reported in that comparison deviate from the values obtained here, with the typical minima 1% lower than the Blacklight value and the typical maxima 1% higher.

It is already common practice to treat any part of a simulation with σ ≳ 1 as vacuum, given that GRMHD jets are often numerically mass loaded and would be unrealistically bright if simulated values there were taken at face value. In fact, this cut is employed in all simulation images shown in this work. Given the importance of probing the origin of image features, Blacklight offers a simple user interface for making similar cuts on a variety of local conditions. Upper and lower cutoffs can be specified for ρ, n_e, p, Θ_e, B, σ, and/or plasma β⁻¹.

As an illustration, the top panels of Figure 5 show the fiducial image and the same image made by considering β⁻¹ < 1/2 to be vacuum, highlighting only the coronal regions of the simulation. Some parts of the image become brighter with this cut, indicating the full data set results in some foreground absorption by less magnetized plasma.

Zoom In Zoom Out Reset image size

Blacklight also provides users with the ability to cut the data based on position. Users are provided with upper and lower cutoffs on radial coordinate r. They can also cut data on one side or another of an arbitrary plane. The plane is defined to be all points x such that $({x}^{a}-{x}_{0}^{a}){n}_{a}=0$ for a user-provided origin x₀ and normal n. The lower panels of Figure 5 show what happens when this plane is chosen to bisect the domain into a foreground and a background relative to the camera.

While having gravity bend geodesics with respect to even a Cartesian coordinate system is critical for obtaining accurate images of light emitted by the modeled system, the resulting mapping from the 3D space of a simulation to a 2D image is very non-intuitive. As an analysis tool intended to be more general than a mere imager, Blacklight allows the user to toggle the effect of gravity. The resulting images more closely correspond to the projections of data one might obtain with standard 3D visualization software, showing features "where they are" rather than where their light reaches the camera.

There is ambiguity in how to treat spacetime as flat in this way. Blacklight uses the following procedure. Simulation data are transformed into Cartesian Kerr–Schild coordinates, both in terms of fluid elements' positions and in terms of the vector components. That is, we take the simulation to provide, perhaps in a roundabout way, ρ, p or κ_e, u^α , and b^α , all as functions of x^α . These values are then reinterpreted to be functions of Cartesian Minkowski coordinates, with vector components interpreted in that coordinate system. Geodesics are calculated in flat spacetime, and the radiative transfer equation integrates the reinterpreted quantities in flat spacetime as well. Qualitatively similar results could be obtained by, e.g., transforming the simulation to spherical Kerr–Schild and then reinterpreting values as being in spherical Minkowski coordinates.

The effect on images of artificially assuming a flat spacetime is illustrated in Figure 6. Though this technique does not produce images as would be observed in nature, it can be useful to identify which image features arise as a result of lensing. In the example shown, one can clearly see that the black hole shadow is magnified by gravity. This feature is particularly useful when used in conjunction with false-color renderings (Section 4.7 below).

Zoom In Zoom Out Reset image size

Most ray tracing is calculated with the fast-light approximation, whereby each time slice from a simulation is taken, in turn, to represent the system in a time-invariant state while light propagates from the plasma to the camera. This is sufficient for many purposes, but the emitting region can be tens or hundreds of gravitational times across, while the dynamical times can be less than 10 gravitational times, with turbulence and caustic-crossing timescales even shorter. Thus certain applications, especially those concerned with high-frequency variability, require slow-light ray-tracing algorithms that account for the time evolution of the system.

Blacklight can read a sequence of simulation outputs, sampling the entire spacetime volume needed to make an image. Whenever a ray needs simulation values from a given spacetime point, those values are interpolated in time just as in space. Users have the option to use either linear interpolation between the two adjacent time-slices, or else simple nearest-neighbor sampling. The latter can be useful for simulation dumps that are not too finely sampled (Δt ≳ 1 GM/c³), as linear interpolation on timescales longer than the magnetic field coherence time tends to systematically underestimate the magnitude of the field (Dexter et al. 2010).

Figure 7 shows the intensity from a fast-light calculation, together with a comparable slow-light image. The fast-light case is the fiducial simulation at time 10,500 GM/c³, though the camera center (and thus outer cutoff for data being used) is moved in to r = 25 GM/c². Even with such a small spatial extent, the longest-duration geodesics need 159 GM/c³ to propagate from their termination points to the camera. The slow-light calculation uses 160 data dumps, equally spaced from 10,371 to 10,530 GM/c³. The offset in this latter time (the time at the camera) from the fast-light time roughly accounts for the light-travel time for the bulk of the emission to the camera. The two images are quite similar, though the approaching plasma on the left side of the image is able to emit light along a number of rays at different times as it moves relativistically.

Zoom In Zoom Out Reset image size

It is sometimes desirable to use a large number of rays (i.e., pixels) in one region of an image, resolving small-scale features, while not wasting resources in other parts of the image. More codes are adopting adaptive ray-tracing techniques in order to facilitate more efficient computation; Wong (2021) refines regions with particularly long geodesic path lengths, Gelles et al. (2021) use an estimate of the error derived from the intensities of nearby pixels, and Davelaar & Haiman (2022) use a scheme very similar to what we describe here in order to trace through a binary black hole system.

The adaptivity in Blacklight is reminiscent of mesh refinement in Athena++. Rather than considering a pixel to be a point, the code considers a pixel to be a square region of the image plane, sampled by a ray placed at its center. This is illustrated in Figure 8. It is a simple matter to recursively subdivide pixels into four, with each iteration augmenting a given ray with four new rays. The root grid, which we will number adaptive level 0, is domain-decomposed into an integer number of square blocks of pixels, whose size is specified by the user. Each block is analyzed as a whole, though with no information from other blocks, as the code runs. If a refinement condition is triggered, every pixel in the block is divided into four, and the four new blocks covering the region, considered to be at the next adaptive level, are ray-traced and evaluated for further refinement. This process continues up to a user-specified maximum refinement level, and there are no limits to how rapidly refinement can change across the image plane. The code supports single-pixel blocks, allowing fine-grained control over which regions are refined.

Zoom In Zoom Out Reset image size

Blacklight supports two types of refinement criteria. One, similar to static mesh refinement in Athena++ and other hydrodynamical codes, specifies regions in the image plane and corresponding minimum refinement levels. Any block whose center lies within a region and whose level is less than the given minimum will be refined. The other class of criteria is based on evaluating the set of intensity values ${I}_{\nu }^{i,j}$ (i and j indexing pixel column and row) at a chosen frequency across all pixels in the block. Users can specify thresholds in intensity ${I}_{\nu }^{i,j}$ itself; absolute pixel-to-pixel intensity gradient

$$\begin{eqnarray}&&{\left|{{\rm{\nabla }}}_{\mathrm{abs}}{I}_{\nu }\right|}^{i,j}={\left({{\rm{\Delta }}}_{i}^{2}+{{\rm{\Delta }}}_{j}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 38a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}=\displaystyle \frac{1}{2}({I}_{\nu }^{i+1,j}-{I}_{\nu }^{i-1,j}),\end{eqnarray} \tag{ 38b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{j}=\displaystyle \frac{1}{2}({I}_{\nu }^{i,j+1}-{I}_{\nu }^{i,j-1});\end{eqnarray} \tag{ 38c }$$

relative pixel-to-pixel intensity gradient

$$\begin{eqnarray}&&{\left|{{\rm{\nabla }}}_{\mathrm{rel}}{I}_{\nu }\right|}^{i,j}={\left({{\rm{\Delta }}}_{i}^{2}+{{\rm{\Delta }}}_{j}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 39a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}=\displaystyle \frac{2({I}_{\nu }^{i+1,j}-{I}_{\nu }^{i-1,j})}{{I}_{\nu }^{i-1,j}+2{I}_{\nu }^{i,j}+{I}_{\nu }^{i+1,j}},\end{eqnarray} \tag{ 39b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{j}=\displaystyle \frac{2({I}_{\nu }^{i,j+1}-{I}_{\nu }^{i,j-1})}{{I}_{\nu }^{i,j-1}+2{I}_{\nu }^{i,j}+{I}_{\nu }^{i,j+1}};\end{eqnarray} \tag{ 39c }$$

absolute pixel-to-pixel Laplacian

$$\begin{eqnarray}&&{\left|{{\rm{\nabla }}}_{\mathrm{abs}}^{2}{I}_{\nu }\right|}^{i,j}=\left|{{\rm{\Delta }}}_{i}+{{\rm{\Delta }}}_{j}\right|,\end{eqnarray} \tag{ 40a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}={I}_{\nu }^{i-1,j}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i+1,j},\end{eqnarray} \tag{ 40b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{j}={I}_{\nu }^{i,j-1}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i,j+1};\end{eqnarray} \tag{ 40c }$$

and/or relative pixel-to-pixel Laplacian

$$\begin{eqnarray}&&{\left|{{\rm{\nabla }}}_{\mathrm{rel}}^{2}{I}_{\nu }\right|}^{i,j}=\left|{{\rm{\Delta }}}_{i}+{{\rm{\Delta }}}_{j}\right|,\end{eqnarray} \tag{ 41a }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}=\displaystyle \frac{4({I}_{\nu }^{i-1,j}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i+1,j})}{{I}_{\nu }^{i-1,j}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i+1,j}},\end{eqnarray} \tag{ 41b }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{j}=\displaystyle \frac{4({I}_{\nu }^{i,j-1}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i,j+1})}{{I}_{\nu }^{i,j-1}-2{I}_{\nu }^{i,j}+{I}_{\nu }^{i,j+1}}.\end{eqnarray} \tag{ 41c }$$

One-sided differences are substituted as appropriate near the edges of blocks. If a user-specified critical fraction of pixels in a block surpasses a given threshold, the block is refined.

As an illustration, we consider the fiducial simulation with a root resolution of 32² pixels, broken into blocks of 8² pixels. The top panels of Figure 9 show the intensity, optical depth, and circular polarization at root level. Though rapidly computed, these images fail to show any features beyond the presence of a crescent with a spur. We then add up to three levels of refinement, refining any blocks for which $\left|{{\rm{\nabla }}}_{\mathrm{abs}}{I}_{\nu }\right|\gt 5\times {10}^{-5}\ \mathrm{erg}\,{{\rm{s}}}^{-1}\ {\mathrm{cm}}^{-2}\ {\mathrm{sr}}^{-1}\ {\mathrm{Hz}}^{-1}$ for five or more pixels. As shown in the lower panels, this criterion focuses pixels in the bright parts of the image with plentiful features, keeping resolution low in the extremities of the image and in the darker parts of the crescent.

Zoom In Zoom Out Reset image size

As discussed in Section 2.6, Blacklight supports thermal, power-law, and kappa distributions for the electrons. A majority of literature results to date assume purely thermal electrons, though this is largely a de facto choice that should become less common as small-scale plasma research continues to inform the GRMHD community about the properties we expect the electrons to have. As more observations are conducted, especially outside the millimeter, alternative distributions may become more necessary to match what is seen in nature.

Figure 10 shows the fiducial image, made with a thermal distribution whose temperature varies across the simulation, together with nonthermal images. The power law is taken to have p = 3, ${\gamma }_{\min }=4$, and ${\gamma }_{\max }=1000;$ we parameterize the kappa distribution with κ = 4 and w = 1/2. These parameters result in total fluxes comparable to the thermal case, though the images are more extended and diffuse. The lower-right panel shows the results from declaring that the electrons are equally partitioned among the three cases, achieved by simply setting the relative fractions in an input file. While not physically motivated, this demonstrates how the code can be used to quickly investigate the effects of more complicated scenarios, such as adding some power-law electrons to a mostly thermal population.

Zoom In Zoom Out Reset image size

One of the primary goals of Blacklight is facilitating science by connecting the physical processes in a simulation with the resulting observables. One way this is accomplished is by providing a number of auxiliary "images" of quantities other than Stokes parameters. This allows for quick access to properties of geodesics and representative simulation values along those geodesics as functions of image-plane location.

One pair of simple parameters is the time and length of each geodesic. The time here is the difference in Kerr–Schild t between the origin of the geodesic (its termination point when tracing backward) and the camera. For length, Blacklight reports the Kerr–Schild proper length s obtained by integrating Equation (21) along the geodesic. Figure 11 shows these two quantities for the fiducial simulation. The time is particularly useful for gauging how many snapshots are required to perform a slow-light calculation with a given camera.

Zoom In Zoom Out Reset image size

We consider three additional quantities indicative of some total property of each geodesic. The code can provide images of the elapsed affine parameter λ, scaled arbitrarily as defined in Section 2.1; it can report the integrated unpolarized emissivity

$$\begin{eqnarray}&&{I}_{j}=\int {j}_{I}\,{d}\lambda ;\end{eqnarray} \tag{ 42 }$$

and it can calculate the total unpolarized absorption optical depth,

$$\begin{eqnarray}&&\tau =\int {\alpha }_{I}\,{d}\lambda .\end{eqnarray} \tag{ 43 }$$

These three quantities are displayed in the top panels of Figure 12.

Zoom In Zoom Out Reset image size

More useful are averages of simulation quantities motivated by these three variables. As discussed in reference to data cuts (Section 4.1), Blacklight tracks the values of ρ, n_e, p, Θ_e, B, σ, and β⁻¹. For each such quantity q, it can report an affine-parameter-averaged value

$$\begin{eqnarray}&&\langle q{\rangle }_{\lambda }=\displaystyle \frac{1}{\lambda }\int q\,{\rm{d}}\lambda ,\end{eqnarray} \tag{ 44 }$$

which weights the quantity relatively uniformly over the geodesic. Alternatively, the emission-weighted average

$$\begin{eqnarray}&&\langle q{\rangle }_{{I}_{j}}=\displaystyle \frac{1}{{I}_{j}}\int {{qj}}_{I}\,{d}\lambda \end{eqnarray} \tag{ 45 }$$

emphasizes those parts of the ray in which the most light is being emitted. Quantities similar to this have already been put to use in the literature. However, emission weighting can be deceptive when a bright emitting region is masked by a cooler, optically thick part of the simulation, in which case the average is not representative of the conditions where the observed photons are emitted. For this reason, we introduce an optical-depth-integrated value 〈q〉_τ , calculated as the solution to

$$\begin{eqnarray}&&\displaystyle \frac{{d}\langle q{\rangle }_{\tau }}{{d}\tau }=q-\langle q{\rangle }_{\tau }\end{eqnarray} \tag{ 46 }$$

in the same way that I is the solution to the unpolarized radiative transfer Equation (25). As the ray is integrated toward the camera, 〈q〉_τ is pulled toward the local value of q at a rate proportional to the absorptivity α_I . The values of 〈Θ_e〉_λ , $\langle {{\rm{\Theta }}}_{{\rm{e}}}{\rangle }_{{I}_{j}}$, and 〈Θ_e〉_τ are shown in the lower panels of Figure 12.

The auxiliary images of Section 4.6 all correspond to a single scalar value generated for each pixel in the image, not unlike true intensity images. Blacklight supports another type of image, which we call a false-color rendering, in which each ray has a human-perceivable color that is modified as it propagates through the simulation.

Consider a quantity q₁, drawn from the same set available to data cuts (Section 4.1) and auxiliary images (Section 4.6): {ρ, n_e, p, Θ_e, B, σ, β⁻¹}. Associate to this quantity a color c₁; an optical depth per unit length τ_s,1; and threshold values q_1,− and q_1,+, either of which may be infinite. A ray, initialized to black, can have its color c be modified in much the same way as the optical-depth integration of Equation (46),

$$\begin{eqnarray}&&\displaystyle \frac{{d}c}{{d}{\tau }_{1}}={c}_{1}-c,\end{eqnarray} \tag{ 47 }$$

where here the appropriate optical depth is related to proper distance s via

$$\begin{eqnarray}{\tau }_{1}=\left\{\begin{array}{ll}s{\tau }_{s,1}, & {q}_{1,-}\lt q\lt {q}_{1,+};\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 48 }$$

In this way, the ray's color approaches c₁ whenever it passes through regions of spacetime where q is between its cutoff values.

Alternatively, let quantity q₂ have color c₂, opacity α₂, and threshold q_2,0. Then whenever the ray crosses a surface q₂ = q_2,0 (possibly only in a user-specified direction), we can blend its color with c₂ according to

$$\begin{eqnarray}&&c\leftarrow (1-{\alpha }_{2})c+{\alpha }_{2}{c}_{2}.\end{eqnarray} \tag{ 49 }$$

Renderings defined in this way can reveal three-dimensional structure in simulation data. Note that these are more akin to volume renderings in scientific visualization software than to externally illuminated isosurface renderings.

We pause to note that it is all too common to find color blending schemes that work in RGB color space. However, the mixing of colored light being emulated by the above equations will fail to match human color perception if done componentwise in this space, obfuscating the scientific accuracy of the final image. Internally, Blacklight uses the correct CIE XYZ color space designed for this purpose (Smith & Guild 1932), which has the added benefit of having a gamut that includes all colors visible to the human eye. Only the final combined image, generated from an arbitrary number of quantities q_i and their associated parameters, need be projected into RGB space. ⁸

Figure 13 shows an example of a three-component rendering from the fiducial simulation. The camera is placed close to the midplane (θ = 70°), and the field of view is enlarged for the first set of images. Blue traces regions with ρ > 5 × 10⁻¹⁸ g cm⁻³, with ${\tau }_{s}={(30\ {GM}/{c}^{2})}^{-1};$ red traces regions with σ > 1 (the standard σ cut has been removed), with ${\tau }_{s}\,={(25\ {GM}/{c}^{2})}^{-1};$ and green tracks β⁻¹ = 1 surfaces, with α = 0.25. Material outside r = 45 GM/c² is ignored here.

Zoom In Zoom Out Reset image size

The first panel, made assuming flat spacetime, is akin to what general-purpose scientific visualization software would produce. This shows the spatial structure of the simulation data in an intuitive way. The second panel more directly connects the simulation variables to lensed images, as it is made with rays that follow the actual null geodesics of the spacetime. The last two panels illustrate the same dichotomy, but zoomed in to the black hole, with material outside r = 20 GM/c² excluded. Additionally, for these latter renderings, we exclude all material with coordinate z > −r_hor. With this cut, we can clearly tell that much of the highly magnetized plasma seen above the black hole actually resides below it.

Figure 14 provides another illustration of what can be learned from false-color renderings. The central panel shows the true 230 GHz image of the fiducial simulation as seen from the pole (θ = 0), with all material outside r = 10 GM/c² omitted. On the left, we illuminate dense (ρ > 2 ×10⁻¹⁷ g cm⁻³) regions, with ${\tau }_{s}={(5\ {GM}/{c}^{2})}^{-1}$. On the right, we illuminate regions with relativistically hot electrons (Θ_e > 20), with ${\tau }_{s}={(1\ {GM}/{c}^{2})}^{-1}$. Comparisons of this sort can show at a glance what simulation variables are most correlated with image structures. Here, the spiral-armed nature of the image outside the photon ring is matched in density but not electron temperature.

Zoom In Zoom Out Reset image size

The vast majority of ray-traced black hole accretion images in the literature are monochromatic, showing I_ν at a single frequency. It is, however, a simple matter to simultaneously calculate images for a number of different frequencies given a fixed camera and fixed simulation data, and there is information contained in the relationships between images at different frequencies.

Here we explore a novel means of presenting such multifrequency data via true-color images, as seen in an appropriately boosted reference frame. The resulting images are particularly useful for outreach and public engagement, where it is often not emphasized that the ubiquitous black–reddish-orange–yellowish-white color maps used in science are merely imparting eye-pleasing colors on what is essentially gray-scale information.

Given values I_ν,i at multiple frequencies ν_i , we can self-consistently Lorentz-boost the intensity into the human visible range. Suppose we wish frequency ν₀ to boost to wavelength ${\lambda }_{0}^{\prime}$ in a new frame. Then the intensities transform as

$$\begin{eqnarray}&&{I}_{\lambda ,i}^{\prime} ={c}^{4}{\left({\nu }_{0}{\lambda }_{0}^{\prime} \right)}^{-5}{\nu }_{i}^{2}{I}_{\nu ,i}.\end{eqnarray} \tag{ 50 }$$

Given I_λ at sample wavelengths, it is straightforward to integrate against tabulated matching functions (Stockman et al. 1999; Sharpe et al. 2005, 2011) based on modern measured cone response functions (Stockman & Sharpe 2000), yielding the XYZ color coordinates of the spectrum. A standard transformation can project the resulting pixel color into RGB values. The only remaining free parameter is the overall normalization of the XYZ values, which changes the intensity but not the hue of the color. It is essentially (the reciprocal of) the camera's exposure.

We apply this procedure to the fiducial simulation, viewed face-on (θ = 0), at an angle (θ = 45°), and edge-on (θ = 90°). In each case, we use 12 frequencies, evenly spaced in wavelength from 152 to 324 GHz. Choosing ν₀ = 230 GHz and ${\lambda }_{0}^{\prime} =550\ \mathrm{nm}$ (near the center of human perception), the boosted wavelengths $\lambda ^{\prime}$ cover the visual range of 390–830 nm. This corresponds to traveling toward the source at a Lorentz factor of approximately 1200. The resulting images are shown in the upper panels of Figure 15, with the lower panels showing the standard 230 GHz monochromatic brightness maps from the same points of view.

Zoom In Zoom Out Reset image size

Such boosted true-color images of synchrotron emission tend to be largely white in color. This is to be expected, as the spectrum is broad compared to the factor of ∼2 range in sensitivity of the human eye, and it tends to be monomodal. Shifting the peak frequencies to wavelengths longer or shorter than ∼550 nm can make it redder or bluer, but it would be surprising to find more exotic colors as a result of this technique. Still, there are subtle color effects in Figure 15, including the fact that the approaching side of the disk is literally blueshifted in hue.