General Purpose Ray Tracing and Polarized Radiative Transfer in General Relativity

The most fundamental object in the Arcmancer library is MetricSpace<n,Signature>, representing an n-dimensional (semi-)Riemannian manifold of a given signature. Defining a new MetricSpace requires the specification of dimensionality n, one or more charts, and functions returning the components of the metric tensor field and its derivatives in each chart. For four-dimensional semi-Riemannian spaces, the metric signature must also be specified. Arcmancer supports both timelike $(+---)$ and spacelike $(-+++)$ signatures.

A chart is represented as a class chart, which in the current implementation only contains a description and serves to give meaning to a tuple of coordinate numbers. The points on the manifold are implemented as a class ManifoldPoint<MetricSpace>. These can be constructed by specifying n coordinates and the corresponding chart. After this, the components of the point can be requested in any available chart, and the object itself behaves much like the mathematical idea of a point on a manifold (see Appendix A.1).

To transform the components of tensorial objects, the transition functions and their Jacobians between the charts must also be specified. For N different charts, this would naively require $N(N-1)$ transition functions and Jacobians to be implemented. The amount of work increases quadratically. However, when the domains of charts ${\phi }_{i}$, ${\phi }_{j}$, and ${\phi }_{k}$ overlap suitably, the transition function from i-coordinates to j-coordinates fulfills

$$\begin{eqnarray}&&{\phi }_{j}\circ {\phi }_{i}^{-1}={\phi }_{j}\circ {\phi }_{k}^{-1}\circ {\phi }_{k}\circ {\phi }_{i}^{-1},\end{eqnarray} \tag{ 1 }$$

and the Jacobian ${J}_{i\to j}=d({\phi }_{j}\circ {\phi }_{i}^{-1})$ decomposes similarly,

$$\begin{eqnarray}&&{J}_{i\to j}={J}_{k\to j}{J}_{i\to k}.\end{eqnarray} \tag{ 2 }$$

The Arcmancer library uses the properties (1) and (2) to build a directed graph of charts, wherein each chart is a node, and the Jacobians and transition functions define the edges. This makes it possible to introduce N charts while supplying only the minimum number of $2(N-1)$ transition functions and Jacobians to make the graph connected. Then, when the components of a point or a tensorial quantity are requested in a different chart, the code walks through the graph, building the transition function and Jacobian piece by piece using Equations (1) and (2).

For a listing of the built-in metric spaces and chart implementations provided with Arcmancer, see Appendix C.

The tensor algebra and calculus for tensors of arbitrary rank (see Appendix A.2) are provided by the Tensor<MetricSpace,Indices...> template class. Here, MetricSpace is the base manifold and Indices... is an arbitrary combination of index tags Cov and Cnt, for covariant or contravariant index, respectively. The implementation is pointwise, using a set of ${n}^{k+l}$ components in a given chart to specify a tensor of rank (k, l) on a manifold M with $\dim (M)=n$ at a given ManifoldPoint.

As such, similarly to a ManifoldPoint, defining a tensor at a point requires the input of ${n}^{k+l}$ components and the corresponding chart. After this, the chart is abstracted away in the sense that algebraic operations between tensors defined at the same point can be performed irrespective of the chart the tensors were originally defined in. The Tensor class provides all of the usual algebraic tensor operations: sum of tensors of the same rank, tensor product, contraction, and in addition, raising and lowering of the indices using the underlying MetricSpace structure. The implementation checks all operations for index correctness at compile time so that, e.g., no contraction between indexes of same type is allowed. In addition, during runtime, all operands are inspected to ensure that they are defined at the same base point. These checks guarantee that operations expressed in code correspond to mathematical operations that are well defined.

The Tensor class also provides some elements of tensor calculus. Namely, the class automatically computes the derivatives required for parallel transporting a tensor along a general curve. Given a curve tangent vector u^a, the class can compute the contractions with ${{\rm{\Gamma }}}_{{bc}}^{a}{u}^{c}$ required in the parallel transport Equation (39).

Functionality for working with curves γ, including geodesics, is provided by the class ParametrizedCurve<MetricSpace,TransportedType> along with a convenience subclass Geodesic. Curves are implemented as sequences of points $(\lambda ,p)$ on a manifold, where λ is the curve parameter and $p\in M$. More concretely, the implementation is based on an ordered queue of objects of type ParametrizedPoint<MetricSpace,TransportedType>, which combines a ManifoldPoint with a real value λ specifying the position along the curve. In addition, the ParametrizedPoint can include any arbitrary object A of type TransportedType to be parallel transported along a geodesic or Fermi–Walker transported along a forced curve. The only requirement is that the object be representable as a (chart-dependent) tuple of real numbers and that a function ${D}_{A}({u}^{a},{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{c},{f}^{a})$ yielding the derivatives $\displaystyle dA(\gamma (\lambda ))/d\lambda$ is provided. The function D_A depends externally on the current tangent vector of the curve u^a, the contractions ${{\rm{\Gamma }}}_{{bc}}^{a}{u}^{c}$, and optionally the force f^a. As mentioned above, the Tensor class provides the derivative function automatically, and as such, arbitrary tensors can be parallel transported along all generic curves without any extra programming effort.

In practice, a curve is computed by specifying the initial conditions in some given chart. These consist of the initial point $({\lambda }_{0},{p}_{0})\in {\mathbb{R}}\times M$, the components of the curve tangent vector ${u}^{a}({p}_{0})\in {T}_{{p}_{0}}M$, the components $A({p}_{0})$ of the possible parallel transported object, and an optional force function f^a. The Arcmancer library then computes points along the curve for the desired interval $I\subset {\mathbb{R}}$ containing ${\lambda }_{0}$ by solving the set of equations (see Appendix A.4)

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma (\lambda )}{d\lambda }={u}^{a},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}^{a}(\lambda )}{d\lambda }=-{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{b}{u}^{c}+{f}^{a},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dA}(\lambda )}{d\lambda }={D}_{A}({u}^{a},{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{c},{f}^{a})\end{eqnarray} \tag{ 5 }$$

in a suitable chart (see Section 3.5 for details on the chart selection).

Arcmancer computes the solution using the integration methods offered by the Odeint C++ library (Ahnert & Mulansky 2011). The default method is the Dormand–Prince fifth-order Runge–Kutta method (Dormand & Prince 1980), which offers error estimation and automatic step-size adjustment, as well as a fair numerical performance in most cases. The absolute and relative error tolerances and step-size and iteration limits are fully user-configurable. After the computation is finished, the ParametrizedCurve class provides access to the solution in any chart and for any $\lambda \in I$. Internally, this is achieved through a cubic spline interpolation.

An interface for implementing hypersurfaces is available through the class Surface<MetricSpace>. Surfaces are useful for representing solid or highly optically thick objects, or regions of interest. Examples include not only the surfaces of neutron stars, white dwarfs, or planets, but also black hole event horizons, optically thick accretion disks, or the limits of computational domains. The Arcmancer implementation of surfaces is based on the concept of level hypersurfaces (see Appendix A.5).

A new surface is implemented by supplying a real-valued function S taking a ManifoldPoint as an argument, as well as the gradient ${\partial }_{a}S$. The surface is then defined as the set of points $\{p\in M| S(p)=0\}$. In addition, a tangent vector field t^a on the surface must be defined. This field is primarily used to represent the four-velocity field of observers fixed on the surface and is required for, e.g., computations involving rotating neutron stars (see Section 6.2).

The Arcmancer library automatically detects intersections of curves with surfaces and numerically finds the exact (to within tolerance) intersection point. The intersections are found by examining the sign of the product $S({p}_{k+1})S({p}_{k})$ for two successive points ${p}_{k+1}$ and p_k on a curve. If the product is negative, the two points must lie in different regions bounded by the surface. The exact intersection point is then found using the so-called Henon's trick (Henon 1982). The "trick" consists of changing the independent variable γ, the curve parameter, in Equations (3)–(4) to S, or the value of the surface function. The transformed equations read

$$\begin{eqnarray}&&\displaystyle \frac{d\lambda }{{dS}}={({u}^{b}{\partial }_{b}S)}^{-1},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma (S)}{{dS}}={({u}^{b}{\partial }_{b}S)}^{-1}\,{u}^{a},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}^{a}(S)}{{dS}}={({u}^{b}{\partial }_{b}S)}^{-1}(-{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{b}{u}^{c}+{f}^{a}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dA}(S)}{{dS}}={({u}^{b}{\partial }_{b}S)}^{-1}\,{D}_{A}({u}^{a},{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{c},{f}^{a}).\end{eqnarray} \tag{ 9 }$$

These equations can then be numerically propagated for a single step of length $-S({p}_{k+1})$ starting from the point ${p}_{k+1}$ to yield the intersection point to within numerical tolerance.

Perhaps the most novel and interesting feature of Arcmancer is the possibility to use multiple coordinate charts simultaneously and seamlessly. The most immediate benefit is that objects can be input and output in any available chart, with all transformations handled automatically by Arcmancer. However, there are important computational benefits to the free selection of coordinate charts as well. The most obvious benefit is the fact that a given problem may be much easier to solve numerically in some specific coordinates compared to others. This is illustrated in Figure 2, where the same null geodesic in an extremal Kerr spacetime is shown in the outgoing Kerr–Schild (KS) coordinates, the ingoing KS coordinates, and the Boyer–Lindquist (BL) coordinates (see Appendix C.2.2). From the figure, it is easy to appreciate how in the outgoing KS coordinates the geodesic is essentially straight, and long integration steps can be taken. On the other hand, in the ingoing KS coordinates and the BL coordinates, the geodesic twists around the event horizon at an increasing rate as the event horizon is approached. The magnitudes of the derivatives with respect to the curve parameter increase correspondingly, making the problem eventually numerically impossible to solve.

Zoom In Zoom Out Reset image size

The possibility to simultaneously use multiple charts makes it possible to avoid the coordinate singularities present in any single chart, such as the pole singularity in any spherical coordinate system or the coordinate singularity at the event horizon present in the usual Schwarzschild coordinates. In addition, using multiple charts makes it possible to switch the chart used for solving the equations of motion for a curve on the fly, useful for situations such as the one depicted in Figure 2. It is not obvious which chart is to be preferred, which is why Arcmancer currently implements several heuristics for automatically choosing the numerically optimal chart.

The first heuristic consists of finding a chart ${\phi }_{i}$ where the matrix of the components of the metric ${\boldsymbol{G}}=({g}_{{ab}})$ has the largest inverse condition number, defined as the ratio of the smallest and largest singular values of the matrix, i.e., ${\mathrm{cond}}^{-1}({\boldsymbol{G}})={\sigma }_{\min }/{\sigma }_{\max }$. This is based on two key observations. First, floating point addition and subtraction between numbers of different magnitude cause a loss of precision. Second, the equations of motion for a curve and for parallel transport along it, Equations (3)–(5), contain a mix of the components of the metric and its derivatives on the right-hand side. As such, it would be intuitively advantageous to perform the computations in a chart where the matrix formed by the metric has eigenvalues that span as small a range as possible. This is achieved by maximizing the inverse condition number.

In some cases, the condition number of the metric is not enough to detect a computationally awkward chart. For example, in the case of a Kerr black hole, the condition number cannot differentiate between the ingoing and outgoing KS charts. However, as is seen in Section 5.1.2, using one over the other can cause a large difference in computation time and accuracy for radial geodesics, depending on whether they are falling toward or emanating from the event horizon. As such, a further heuristic is needed.

If the condition number heuristic does not separate two promising charts, the Arcmancer code next tries to minimize the maximal absolute value of the intrinsic derivatives, $-{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{b}{u}^{c}$, of the curve tangent vector u^a. As such, this heuristic needs to know the current curve tangent vector u^a, unlike the condition number test, for which only the current point is required. For Cartesian coordinates in a Euclidean or Minkowskian space, ${{\rm{\Gamma }}}_{{bc}}^{a}\equiv 0$, so in effect this procedure looks for the most Cartesian-like chart in which the metric looks most Euclidean (or Minkowskian) in the direction of the current curve tangent vector u^a.

Formal proofs of the performance of these heuristics are beyond the scope of this work, but the numerical results in Section 5 indicate that they work reasonably well.

For four-dimensional Lorentzian manifolds, Arcmancer provides a functionality to construct local Lorentz frames (see Appendix A.6) through the class LorentzFrame<MetricSpace>. The user supplies a timelike vector e_t and two spacelike vectors e_z and e_x. From these, a complete Lorentz frame $\{{E}_{t},{E}_{x},{E}_{y},{E}_{z}\}$ is constructed by first normalizing e_t to yield E_t and then orthonormalizing e_z and e_x sequentially. Finally, E_y is defined by the remaining orthogonal direction through ${E}_{y}^{a}=\pm \epsilon {}_{{bcd}}^{a}{E}_{t}^{b}{E}_{z}^{c}{E}_{x}^{d}$, where ${\epsilon }_{{abcd}}$ is the Levi–Civita tensor, with the sign depending on the desired handedness (positive for a right-handed frame).

The LorentzFrame object can be automatically parallel transported along a ParametrizedCurve. In addition, Tensor objects can be constructed from components given with respect to a LorentzFrame. Likewise, the components of any Tensor can be extracted in a given LorentzFrame as well; see Equation (41).

To produce mock observations, an observational instrument must be emulated somehow. For ray-tracing purposes, this usually means using an image plane. The image plane is positioned near the object of interest, and only the rays intersecting the plane orthogonally are considered. These rays are then assumed to propagate in vacuum all the way to the distant observer. This approximation ignores atmospheric and instrumental effects, but these can be modeled afterwards using dedicated tools if necessary.

There are three main sources of error when generating the image plane: required deviations from perpendicularity, perturbations caused by the curvature of the space, and the assumption of vacuum propagation. The assumption of perpendicularity is typically excellent. For distant objects, the maximum deviation from perpendicularity ${\rm{\Delta }}\theta$ is approximately equal to the observed angular size of the object, or ${\rm{\Delta }}\theta \sim L/(2D)$, where L is the linear extent of the source perpendicular to the line of sight and D is the distance. For example, in the case of Sgr A* (Sagittarius A^*), we have ${\rm{\Delta }}\theta \sim {10}^{-11}$, and for a typical galactic neutron star ${\rm{\Delta }}\theta \sim {10}^{-17}$. The effects of the remaining spacetime curvature at the image plane location can be estimated by looking at the bending angle β that the image plane rays will make when propagated to infinity. Sufficiently far away from the object so that the Schwarzschild metric can be used, this angle turns out to be (e.g., Beloborodov 2002) $\beta \sim 2{GM}/({c}^{2}R)$, where M is the total mass of the observed object and R is the radial distance of the image plane from the object. Thus, for $R\gtrsim {10}^{4}\,{GM}/{c}^{2}$, we have $\beta \lesssim 2\times {10}^{-4}$, and so the effects of residual curvature are negligible. The assumption of propagation in vacuum is typically valid for objects that are not situated at cosmological distances as far as light bending is concerned. However, corrections for effects such as extinction, frequency dispersion, or Faraday rotation may need to be added in further postprocessing.

In many ray-tracing codes, the construction of image planes is achieved by a assuming a flat space and explicitly constructing the starting points and tangent vectors for a planar configuration of geodesics (Broderick 2004; Cadeau et al. 2007; Dexter & Agol 2009; Vincent et al. 2011; Dexter 2016; Chan et al. 2017). Arcmancer provides a general purpose tool for constructing plane-parallel initial conditions for Lorentzian spacetimes in the class ImagePlane<MetricSpace,DataType>. The user specifies a LorentzFrame at the center of the plane and the extent and the resolution (number of grid points) of the plane in the local E_x and E_y directions. The local Lorentz frame is then parallel transported to the desired grid points via spacelike geodesics, using the Arcmancer curve propagation functionality. Initial conditions for curves passing through the plane are set up by assigning the tangent vectors ${u}^{a}(0)$ to be spatially parallel to the parallel transported E_z vector. The collection of parallel transported frames defines a best local approximation to a flat plane that is threaded by orthogonal curves and corrects the effect of the bending β caused by the curvature to first order. Thus, the Arcmancer ImagePlane can safely be used in regions where the curvature is small but non-negligible. The method is also general purpose in the sense that it works similarly in any coordinate system and only requires specifying a local Lorentz frame at one point.

Radiative transfer functionality in Arcmancer is built with flexibility in mind. For this purpose, the interface declares two types of functions. The first type of function is FluidFunction<MetricSpace,FluidData>, which maps points on the base manifold MetricSpace to a user-defined set of fluid variables FluidData, which represent local material properties such as temperature or density. The only restriction is that FluidData must include a single bulk fluid four-velocity w^a and a single reference direction (often magnetic field) t^a orthogonal to w^a.

The second type of function is RadiationFunction<FluidData>, which computes the Stokes emissivity vector ${\boldsymbol{J}}$ and the response matrix ${\boldsymbol{M}}$ (see Appendix B) from the given FluidData, local fluid rest-frame frequency ν, and the rest-frame angle θ between the reference direction t^a and the current direction of the light ray (the tangent vector k^a).

This approach makes implementing different fluid and radiation models rather straightforward. For example, the fluid variables for a given point can be obtained from a GR magnetohydrodynamics (GRMHD) simulation or from an analytic model. The Arcmancer suite includes an example application which reads outputs from the HARM GRMHD code (Gammie et al. 2003; Noble et al. 2006) and computes mock observations using a thermal synchrotron radiation model based on the results in Dexter (2016). See Section 5.2.2 for computational results.

With Arcmancer, a radiative transfer problem (see Appendix B) is solved by first propagating a set of curves ${\gamma }_{i}$ (typically geodesics, unless plasma effects are significant) along which the radiative transfer equation, Equation (52), is to be solved as a curve integral. Usually, the most convenient approach is to use an ImagePlane and let Arcmancer propagate the set of initial conditions backwards in time through the region of interest. Each propagated curve must include a parallel transported PolarizationFrame, a pair of two orthogonal spacelike vectors ${ \mathcal P }=\{{v}^{a},{h}^{a}\}$, also orthogonal to the geodesic and the four-velocity of the observer, representing the vertical and horizontal linear polarization basis vectors of the observer at one end ${\gamma }_{i}({\lambda }_{\mathrm{obs}})={p}_{\mathrm{obs}}$ of the curve. If using an ImagePlane, these can be conveniently obtained from the E_x and E_y vectors of the local Lorentz frame at each point.

The four-velocity ${u}^{a}({p}_{\mathrm{obs}})$ of the observer ${ \mathcal O }$ at ${p}_{\mathrm{obs}}$, the four-velocity ${w}^{a}(p)$ of the fluid at each point $p={\gamma }_{i}(\lambda )$, and the curve tangents ${k}^{a}(p)$ and ${k}^{a}({p}_{\mathrm{obs}})$ define a connection between the photon frequency ${\nu }_{0}$ observed by ${ \mathcal O }$ at ${p}_{\mathrm{obs}}$ and the corresponding photon frequency ν in the local rest frame of the fluid at p. This is given by the redshift factor

$$\begin{eqnarray}&&{ \mathcal G }=\displaystyle \frac{{u}_{a}({p}_{\mathrm{obs}}){k}^{a}({p}_{\mathrm{obs}})}{{w}_{a}(p){k}^{a}(p)}=\displaystyle \frac{{\nu }_{0}}{\nu }.\end{eqnarray} \tag{ 10 }$$

The initial conditions are set by defining initial invariant specific intensities $\{{{\boldsymbol{ \mathcal I }}}_{{\nu }_{0,i}/{ \mathcal G }}={{\boldsymbol{I}}}_{{\nu }_{0}/{ \mathcal G }}{({\nu }_{0}/{ \mathcal G })}^{-3}\}$ at the other end ${p}_{\mathrm{start}}$ of the curve, one for each observed frequency ${\nu }_{0,i}$ of interest. Often these can be set to zero, but, for example, in the case of radiation emanating from optically thick or solid surfaces, the initial intensity can be non-zero. Solving the radiative transfer equation itself proceeds in a manner following Shcherbakov & Huang (2011). See Figure 3 for a diagram of all the vectors and angles.

Zoom In Zoom Out Reset image size

At each point $p\in M$ during the calculation, the Arcmancer library evaluates the given FluidFunction to obtain the fluid four-velocity w^a and the rest of the fluid parameters in the rest frame of the fluid. This includes the local reference direction b^a, which typically is the direction of the local magnetic field. From these, the angle $\theta (w;b,k)$ between the reference direction b^a and the light ray tangent k^a as seen in the fluid rest frame is computed using Equation (45). This angle is required by some radiation models, such as synchrotron emission models. The reference direction also defines the local vertical direction of polarization ${\boldsymbol{V}}={\boldsymbol{k}}\times ({\boldsymbol{k}}\times {\boldsymbol{b}})$, where ${\boldsymbol{k}}$ and ${\boldsymbol{b}}$ are the spatial parts of k^a and b^a, respectively.

The next step is to project the parallel transported polarization frame ${ \mathcal P }$ to the fluid rest frame using the screen projection operator, Equation (44), yielding $\tilde{{ \mathcal P }}=\{{\tilde{v}}^{a},{\tilde{h}}^{a}\}$, where

$$\begin{eqnarray}&&{\tilde{v}}^{a}={P}_{\perp }{(w,k)}_{b}^{a}{v}^{b},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\tilde{h}}^{a}={P}_{\perp }{(w,k)}_{b}^{a}{h}^{b}.\end{eqnarray} \tag{ 12 }$$

Now we can compute the angle χ between the projected parallel transported polarization frame $\{{\tilde{v}}^{a},{\tilde{h}}^{a}\}$ and the polarization frame of the fluid, defined by ${\boldsymbol{V}}$, from

$$\begin{eqnarray}&&\tan \chi =\displaystyle \frac{-{V}_{a}{\tilde{h}}^{a}}{{V}_{a}{\tilde{v}}^{a}},\end{eqnarray} \tag{ 13 }$$

where ${V}^{a}=(0,{\boldsymbol{V}})$.

Next, the angle θ and the fluid parameters are passed to the RadiationFunction to obtain the Stokes emissivity and the response (Müller) matrix ${{\boldsymbol{M}}}_{\nu }$ in the fluid rest frame. These are related to the parallel transported and projected polarization frame $\tilde{{ \mathcal P }}$ using the angle χ and the transformation properties of the Stokes components under rotation (e.g., Chandrasekhar 1960). The emissivity vector ${{\boldsymbol{J}}}_{\nu }$ and response matrix ${{\boldsymbol{M}}}_{\nu }$ are transformed via ${{\boldsymbol{J}}}_{\nu }\mapsto R(\chi ){{\boldsymbol{J}}}_{\nu }$ and ${{\boldsymbol{M}}}_{\nu }\mapsto R(\chi ){\boldsymbol{M}}R(-\chi )$, where

$$\begin{eqnarray}R(\chi )=\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & \cos (2\chi ) & -\sin (2\chi ) & 0\\ 0 & \sin (2\chi ) & \cos (2\chi ) & 0\\ 0 & 0 & 0 & 0\end{array}\right)\end{eqnarray} \tag{ 14 }$$

gives the transformation of Stokes vectors under rotations of the polarization plane. Finally, it can be shown that the Stokes components in any two polarization frames ${ \mathcal P }$ and $\tilde{{ \mathcal P }}$ related by a screen projection are equal, so that the radiative transfer equation to be solved along the geodesic is

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{ \mathcal I }}}_{\nu }}{d\lambda }=L\displaystyle \frac{{u}_{a}({p}_{\mathrm{obs}}){k}^{a}({p}_{\mathrm{obs}})}{{\nu }_{0}}({{\boldsymbol{ \mathcal J }}}_{\nu }-{{\boldsymbol{ \mathcal I }}}_{\nu }{{\boldsymbol{ \mathcal M }}}_{\nu }),\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal J }}}_{\nu }={\nu }^{-2}R(\chi ){{\boldsymbol{J}}}_{\nu },\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal M }}}_{\nu }=\nu R(\chi ){{\boldsymbol{M}}}_{\nu }R(-\chi ),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\nu ={\nu }_{0}/{ \mathcal G },\end{eqnarray} \tag{ 18 }$$

and L is the unit of length. For example, in problems related to black holes, a typical choice is $L={GM}/{c}^{2}$, where M is the black hole mass. Internally, Equation (15) is solved using the Odeint Runge–Kutta–Fehlberg eighth-order method. However, for problems where the optical thickness is large, Equation (15) can become stiff, and an implicit method would provide better performance.

5.1.1. Geodesic Propagation

The accuracy of the basic curve propagation functionality (Section 3.3) was verified by investigating curves on a two-dimensional spherical surface. The computations were performed both in two dimensions, using the intrinsic spherical coordinate chart $(\theta ,\phi )$, Equation (53), and in a three-dimensional Euclidean slice at t = 0 of the Minkowski space, using the spherical coordinates $(0,r,\theta ,\phi )$, Equation (55). To force the curve to stay on the surface of a sphere in the three-dimensional case, a constraint force $f({u}^{a})=(0,-{u}^{a}{u}_{a},0,0)$ was specified. Here, u^a is the curve tangent, in three-dimensional spherical coordinates.

Numerical convergence was estimated using a single geodesic curve $\gamma (\lambda )$ passing through $(\theta ,\phi )=(\pi /2,0)$ at $\lambda =0$ with a tangent vector ${u}^{a}=(\dot{\theta },\dot{\phi })=(1.0,0.3)$. The initial values were chosen so as to avoid a purely polar or equatorial geodesic, but were otherwise chosen arbitrarily. The geodesic was computed several times using a range of equal relative and absolute numerical tolerances ${\epsilon }_{\mathrm{rel}}$ and ${\epsilon }_{\mathrm{abs}}$ from 10⁻²⁰ to 10⁻² in 40 steps. The differences between the numerical results and the known analytical solution are shown in Figure 4. We see that both in the intrinsic two-dimensional and the constrained three-dimensional case, the numerical curves converge toward the analytical solution linearly with the tolerance parameters. The convergence saturates at tolerance parameters $\sim {10}^{-15}$ when the relative precision floor of the double precision floating point numbers is reached.

Zoom In Zoom Out Reset image size

5.1.2. Parallel Propagation in the Kerr Spacetime

The functionality for parallel transporting tensorial quantities (see Sections 3.2 and 3.3) along a curve was assessed in the context of a Kerr spacetime (see Section C.2.2) with a near-extremal non-dimensional spin parameter $\chi =0.99$ and mass M = 1. First, the initial conditions $\gamma (0)=({t}_{0},{r}_{0},{\theta }_{0},{\phi }_{0})$ $=(0,10,1,0)$ and ${u}^{a}(0)=(1,0,0.01,0.03)$ were fixed in the BL coordinates (see Equation (56)). These initial values were chosen to yield a generic timelike geodesic, and to avoid special cases such as equatorial geodesics, but were otherwise chosen arbitrarily. The geodesic was then augmented by including the metric g_ab and a Lorentz frame $\{{E}_{t},{E}_{x},{E}_{y},{E}_{z}\}$ as quantities to be parallel transported. The geodesic was then computed until $\lambda =990$ to yield several complete orbits around the black hole, using tolerances ${\epsilon }_{\mathrm{abs}}={\epsilon }_{\mathrm{rel}}={10}^{-10}$. Finally, the parallel transported values were evaluated for accuracy by comparing to analytic expectations.

Figure 5 shows the orbit of the geodesic. It also depicts magnitudes of the maximum difference $\max | {\rm{\Delta }}{g}_{{ab}}|$ of the components of the parallel transported metric with respect to the analytic expression, both computed in the ingoing KS chart. Also shown are the absolute values of all the pairwise inner products of the parallel propagated Lorentz frame which should be identically zero. From the figure, we see that the errors in all of these conserved quantities increase in a secular fashion, while the single step errors are below the set numerical tolerance. This is an expected and well-known behavior for non-symplectic numerical integration methods, such as the fifth-order Dormand–Prince scheme used in Arcmancer, which does not respect the geometric structure of the phase space (Hairer et al. 2008). Symplectic methods for the inseparable Hamiltonians occurring in geodesic propagation have been discovered recently (Pihajoki 2015), but these are not yet available in Odeint. In general, the secular accumulation of integration error poses no problem for the applications we demonstrate in this paper. However, for integrations over long periods of time, such as for computing dynamics of massive particles orbiting a black hole, a symplectic method for inseparable Hamiltonians might need to be implemented.

Zoom In Zoom Out Reset image size

The accuracy and performance of both the curve propagation and parallel transport functionality were also assessed as a function of the geodesic and the coordinate chart. To this end, we set up an image plane at $({r}_{0}={10}^{5},{\theta }_{0}=50^\circ )$ in the BL coordinates of a Kerr spacetime with $\chi =0.95$ and M = 1. From the image plane, null geodesics were propagated backwards from $\lambda =0$ to $\lambda =-2{r}_{0}$ or until intersection with a surface slightly outside the event horizon, defined by $r=1.03{r}_{H}$, where r_H is the event horizon radius. This radius was chosen since the computation in the BL and ingoing KS coordinates must be terminated before the event horizon itself (see Figure 2). The geodesics were computed three times, each time fixing the chart (automatic chart selection disabled) to either ingoing KS, outgoing KS, or the BL chart. Standard tolerances of ${\epsilon }_{\mathrm{rel}}={\epsilon }_{\mathrm{abs}}={10}^{-10}$ were used. We then computed the maximal absolute errors in the value of the curve Hamiltonian, $H(x,k)={k}_{a}{k}^{a}\,=\,0$, and the trace ${{g}^{a}}_{a}$ of the parallel transported metric along the geodesics and plotted these on the image plane, in addition to the number of integration steps N. The results are shown in Figure 6.

Zoom In Zoom Out Reset image size

From the figure, it is evident that the outgoing KS coordinates offer significantly better numerical performance than the ingoing KS coordinates or the BL coordinates. This is not surprising, since the outgoing KS chart is adapted to radially outgoing null geodesics. As a consequence, the more radial the geodesic is, the more nearly a straight line it is in the outgoing KS chart. In the figure, this can be seen as the remarkable decrease in the maximal error and the number of computational steps for the geodesics starting near the origin of the image plane (see also Figure 2). On the other hand, the BL coordinates are seen to perform significantly worse. This is related to both the fact that the geodesic "wraps around" the black hole near the event horizon (see Figure 2), but also the fact that the condition number (the ratio of the maximum to minimum singular value) of the matrix of the metric components scales as ${ \mathcal O }({r}^{2})$ (see Section 3.5). In addition, the coordinates are singular at the poles. All of these factors combine to make the BL chart the most numerically disadvantageous of the three. Finally, the ingoing KS chart fares worse than the BL chart for the conservation of the Hamiltonian, but better for the trace of the metric tensor and number of steps taken. This is understandable, since these coordinates are adapted to radially ingoing null geodesics, and outgoing geodesics "wrap around" the black hole near the event horizon twice as fast compared to the BL coordinates. This is partly offset by the fact that the condition number of the metric components is better behaved than for the BL coordinates. The "wraparound" behavior is suppressed near the poles of the black hole, which in the figure can be seen as the slight decrease in the error of the metric trace around the "north" pole of the black hole for the BL and the ingoing KS coordinates.

The accuracy in general is seen to be consistent with the given numerical tolerances. The outgoing KS chart in particular provides excellent accuracy, with results much better than even the set tolerances for nearly radial geodesics. In addition, there is a factor of ∼10 difference in the number of steps taken between the outgoing KS chart and the BL chart, which was also directly reflected in the computational time. The results strongly suggest that the outgoing KS metric should be preferred in all codes computing mock observations using geodesics emanating from the vicinity of a Kerr black hole. Likewise, for studies of radiation scattering from a black hole, ingoing KS coordinates should be used for computing the incoming radiation and outgoing KS coordinates for the scattered, outgoing radiation.

We assessed the accuracy and convergence properties of the Arcmancer radiative transfer functionality by postprocessing a general relativistic magnetohydrodynamics (GRMHD) simulation and by comparing to an existing polarized radiative transfer code grtrans (Dexter 2016). To facilitate an easy comparison, we used the same simulation data used to test grtrans, as the data are conveniently distributed with the grtrans code.⁵ The simulation data used were computed with the GRMHD code HARM (Gammie et al. 2003; Noble et al. 2006) and describe an axisymmetric optically thin accretion flow around a Kerr black hole with a dimensionless angular momentum of χ = 0.9375. The black hole mass M and its accretion rate $\dot{M}$ were set to the grtrans defaults for HARM, $M=4\,\times {10}^{6}\,{M}_{\odot }$ and $\dot{M}=1.57\times {10}^{15}\,{\rm{g}}\,{{\rm{s}}}^{-1}$ $=2.49\times {10}^{-11}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Assuming in addition that the source is at a distance of D = 8 kpc, these values approximate Sgr A* (Dexter 2016), although for this particular data set at the observed frequency of ν = 230 Ghz the total computed flux of $\sim 13\,\mathrm{Jy}$ (see below) is roughly three times too large compared to the observed value of ∼3–4 Jy (e.g., Bower et al. 2015). The radiative model was taken to be relativistic thermal synchrotron radiation, using the updated formulae in Dexter (2016). The electron and proton temperatures in the plasma were assumed equal, and the ideal gas equation of state was assumed, also corresponding to grtrans.

All mock observations were computed using a square image plane with physical dimensions $x,y\in [-L,L]$, where $L\,=13\,M$, in the local Lorentz frame of a stationary observer with ${E}_{t}=(1,0,0,0)$ (see Section 3.6 and Appendix A.6.1). The negative z-axis is pointed toward the origin of the BL coordinates. The image plane was set at a distance $r={10}^{4}\,M$ and an inclination of $\theta =50^\circ$ (in BL coordinates), as in Dexter (2016). Geodesics from the image plane were then computed backwards in time from the image plane and the radiative transfer computed at the observed frequency of $\nu =230\,\mathrm{GHz}$ along these geodesics to form the final image. Numerical tolerances for both the geodesic computation and the radiative transfer computation were set to 10⁻¹⁰. These tolerances guarantee that the accuracy during radiation transfer is governed by the chosen sampling rate ${\rm{\Delta }}\lambda /{\lambda }_{\max }$, where ${\lambda }_{\max }$ is the total range of the affine parameter over which the radiation transfer is computed. This ensures that the convergence and comparison results are not affected by the characteristics of the sampling induced by time-step control.

5.2.1. Flux Convergence

First, we investigated the convergence of the total flux in the Stokes variables ${{\boldsymbol{I}}}_{\nu }=({I}_{\nu },{Q}_{\nu },{U}_{\nu },{V}_{\nu })$ as the maximum step size ${\rm{\Delta }}\lambda$ in the affine parameter and the image size P in pixels per side were varied. For each pixel i, we computed the observed flux

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\nu ,i}={{\boldsymbol{I}}}_{\nu }({x}_{i},{y}_{i})\displaystyle \frac{{\rm{\Delta }}x{\rm{\Delta }}y}{{D}^{2}},\end{eqnarray} \tag{ 19 }$$

where ${\rm{\Delta }}x={\rm{\Delta }}y=2L/P$ is the physical size of each pixel and D is the (non-cosmological) distance of the target. From the pixel-by-pixel fluxes, the total integrated fluxes ${{\boldsymbol{F}}}_{\nu }={\sum }_{i}{{\boldsymbol{F}}}_{\nu ,i}$ were then computed.

The convergence results are shown in Figure 7. The general trend is that the benefits of a smaller step size saturate quickly for smaller-sized images, where the spatial sampling noise dominates. Similarly, increasing the image size is only effective up to the point where the noise from sampling of the small-scale structures starts to dominate. For this particular case, the benefit of increasing the image size beyond P = 316 pixels per side is already marginal. At this size, a $0.1 \%$ convergence is achieved at a maximum relative step size of ${\rm{\Delta }}\lambda /{\lambda }_{\max }=3\times {10}^{-3}$.

Zoom In Zoom Out Reset image size

5.2.2. Comparison to grtrans

In addition to ensuring the consistency and convergence of the Arcmancer results, we performed a comparison to a publicly available radiative transfer code grtrans using the same HARM data set as above. Both codes were used to compute a square image 400 pixels wide, as above. grtrans was configured to take 2000 steps, which according to Dexter (2016) should net a relative accuracy for total flux at the 10⁻³ level. Similarly, Arcmancer was constrained to take steps of at most ${\rm{\Delta }}\lambda /{\lambda }_{\max }=3\times {10}^{-3}$, which should guarantee a relative accuracy of better than 10⁻³ by the convergence results above.

The resulting Stokes intensity maps computed with Arcmancer are shown in Figure 8. In addition, Figure 9 shows the relative differences in the Stokes intensities as computed by Arcmancer versus grtrans. In the left panel of Figure 9, we see that the unpolarized intensity predicted by Arcmancer is consistently higher, and there is a clear difference in the polarized results, especially in the Q and U components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The reason for this discrepancy was traced to two separate numerical issues. First, grtrans uses values for the gravitational constant G and Boltzmann constant k_B that were truncated to three significant figures, while in Arcmancer the CODATA 2014 (Mohr et al. 2016) values are used up to the known experimental precision. Second, grtrans uses an approximation for computing the cylindrical Bessel functions used in the relativistic thermal synchrotron radiation model, Equations (B4) and (B14) in Dexter (2016). The grtrans code, as well as some other codes, such as that in Mościbrodzka & Gammie (2018), use first-order approximations for the cylindrical Bessel functions, but at least in this example case the approximations are not always valid throughout.

If the same physical constants and Bessel function approximations are used in Arcmancer, the agreement both in unpolarized and polarized intensities is excellent, as can be seen in the right panel of Figure 9. The unpolarized and polarized total intensities agree with grtrans to a relative level of $\sim {10}^{-3}$, and the pixel-by-pixel errors are below percent level on average. There are a small number of high difference outliers located either at regions where the absolute intensity values are very small or at the strongly lensed rings of emissivity. The former outliers are caused by numerical noise and the latter mainly by spatial sampling noise, since the small-scale structure in these rings is not resolved while simultaneously the emission is highly boosted, amplifying the differences. However, it can be seen from the results in Figure 7 and the numerical total fluxes tabulated in Table 1 that these pixels make no significant difference in the observed integrated fluxes.

Table 1.  Total Integrated Fluxes as Computed with Arcmancer in the Test Scenario Depicted in Figures 8 and 9, Relative to the Values Obtained with grtrans

I	Q	U	V
Arcmancer
1.0083	0.98675	1.1671	1.0054
With grtrans compatibility
1.0009	1.0029	0.99690	1.0023
With no polarization
1.1264	0	0	0

Download table as:  ASCIITypeset image

For an interesting test case, we also ran the same test scenario with all polarization effects disabled. That is, we set ${j}_{\{Q,U,V\}}={\alpha }_{\{Q,U,V\}}={r}_{\{Q,U,V\}}=0$ in ${\boldsymbol{J}}$ and ${\boldsymbol{M}}$ so that only the unpolarized degrees of freedom were propagated. As shown in Table 1, the resulting total flux is $\sim 13 \%$ higher than in the polarized case. This suggests that creating mock observations of unpolarized flux can be misleading if polarization effects are completely ignored.

Arcmancer makes it easy to compute mock observations of emitting surfaces with different user-defined geometries. Here, this feature is demonstrated through a toy model by computing the changes on the observed spectropolarimetric features caused by varying the opening half-angle β of a geometrically thin but optically thick accretion disk around a Kerr black hole.

Often (e.g., Vincent et al. 2011; Bambi 2012; Psaltis & Johannsen 2012; Dexter 2016), a thin disk is modeled in mock observation simulations as an infinitely thin equatorial plane around the black hole. For α-disk models (Shakura & Sunyaev 1973; Novikov & Thorne 1973), this is in many cases a satisfactory approximation. This is because the maximal angle made by the disk photosphere and the symmetry plane is $\beta \sim \arctan (0.2\dot{m})$, computed in the Schwarzschild coordinates, where $\dot{m}=\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ is the black hole accretion rate in units of the Eddington accretion rate ${\dot{M}}_{\mathrm{Edd}}$. However, for an accretion rate of $\dot{m}=0.3$, this maximal angle is already $\beta \sim 4^\circ$, which can be expected to have observable consequences. This is since the maximal β in the Shakura–Sunyaev solution is found at $r=(27/2)M$, where M is the black hole mass, which is in the bright inner region of the disk.

Instead of the geometry of the α-disk model, which has a photospheric surface profile dependent on the accretion rate, we use a disk defined by a hyperbolic surface in the outgoing KS coordinates,

$$\begin{eqnarray}&&S(t,x,y,z)=\displaystyle \frac{{x}^{2}+{y}^{2}-{a}^{2}}{{r}_{\min }^{2}}-\displaystyle \frac{{z}^{2}}{{r}_{\min }^{2}{\tan }^{2}\beta }-1,\end{eqnarray} \tag{ 20 }$$

where $a=M\chi$ is the normalized angular momentum of the black hole, β is the half-opening angle of the hyperboloid, and ${r}_{\min }$ sets the inner boundary of the disk, here fixed to the innermost stable circular orbit (ISCO) of the black hole. The choice of this surface is motivated by the intention to investigate only the effects of geometry on the observable properties while keeping the emission properties of the disk otherwise fixed.

To compute the mock observation, we first set up an image plane with physical dimensions $x,y\in [-40M,40M]$ at a distance of $r={10}^{4}M$ (in BL coordinates). From this surface, null geodesics were propagated backwards until they intersected the disk surface or the event horizon. A PolarizationFrame was parallel transported with the geodesic to enable Faraday rotation effects to be captured. Points that hit the disk were given a blackbody spectrum ${B}_{\nu }(T)$ with a temperature T matching the Novikov–Thorne disk model (Novikov & Thorne 1973; Page & Thorne 1974), using mass $M=10\,{M}_{\odot }$, accretion rate $\dot{m}=0.3\,{\dot{M}}_{\mathrm{Edd}}$, and dimensionless viscosity parameter $\alpha =0.1$. The intensity and the linear polarization of the point were computed based on the electron scattering atmosphere model given in Chandrasekhar (1960), using the impact angle θ between the geodesic and the disk normal, computed in the rest frame of the rotating disk surface. The exact solution requires solving an integral equation. We instead used the Padé approximants

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\nu }(\mu )}{{S}_{\nu }}=\displaystyle \frac{1+2.3\mu -0.3{\mu }^{2}}{2\times 1.19167}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&P=0.117126\,\displaystyle \frac{1+14.9165\mu -15.8923{\mu }^{2}}{1+22.2420\mu +44.8893{\mu }^{2}},\end{eqnarray} \tag{ 22 }$$

where $\mu =\cos (\theta )$ for the intensity I_ν, normalized by the source function S_ν (in this example, ${S}_{\nu }={B}_{\nu }(T)$), and the polarization fraction P, respectively. Both approximations are accurate to within 2% over the range $\mu \in [0,1]$. It should be noted that the combination of a blackbody spectrum and a beamed intensity profile is not fully self-consistent, since a genuine blackbody emitter is isotropic. However, the combination serves to illustrate the effects of an anisotropically emitting surface. In addition, the spectral shape for thin accretion disks around stellar mass black holes is in any case well described using a diluted blackbody (Davis et al. 2005).

To construct the image from these data, instead of running full radiation transfer, the radiation was assumed to propagate in vacuum. This is not a particularly good assumption physically, since the thin disks are expected to have a tenuous, hot corona (e.g., Liang & Price 1977; Czerny & Elvis 1987), but it was made so as to not add additional uncertainties and keep the focus on the effects of changing disk geometry. The values of the intensity and polarization were directly transferred to the image plane after scaling the intensity by the redshift factor and rotating the polarization to match the rotation of the parallel transported PolarizationFrame. This computation was repeated for seven values of the disk opening angle from β = 0001 to β = 25° and three observer inclination angles i = 10°, 35°, and 60°. The results are collected in Figures 10 and 11, which show mock images of the two extreme cases (β = 0001, β = 25°) and the polarization spectra for all the computation runs.

Zoom In Zoom Out Reset image size

The effect of the disk opening angle is clearly seen in Figure 10, which shows specific intensity maps as seen by an observer at an inclination of $i=60^\circ$ for the extreme opening angles of β = 0001 and β = 25°. The intensity patterns differ significantly, with most of the emission coming from the opposite side of the disk for the disk with the larger opening angle. In addition, the structure of the ring caused by radiation that has traveled around the black hole once is noticeably changed by the increased disk thickness (cf. Luminet 1979). Despite the visual differences, Figure 11 shows that the shape of the spectra obtained from the integrated emission is hardly changed at all, and as such, the shape of the observed spectrum is not very sensitive to the disk geometry in this example.

Zoom In Zoom Out Reset image size

Figure 10 also shows a significant difference in polarization patterns, with the large opening angle disk exhibiting a large asymmetry between the upper and lower halves of the mock observation image. This is caused by a purely geometrical effect, wherein the geodesics emanating from the opposite side of the disk from the observer's point of view are more closely aligned with the local disk surface normal. For the geodesics coming from the observer's side of the disk, the situation is the opposite. The polarization fraction of the electron scattering atmosphere model is strongly dependent on the angle of the geodesic with respect to the disk normal, with stronger polarization for lower incidence angles. The graphs of the degree of polarization,

$$\begin{eqnarray}&&P=\displaystyle \frac{\sqrt{{Q}^{2}+{U}^{2}+{V}^{2}}}{I},\end{eqnarray} \tag{ 23 }$$

and the polarization angle,

$$\begin{eqnarray}&&\psi =\displaystyle \frac{1}{2}\arctan (U/Q),\end{eqnarray} \tag{ 24 }$$

in Figure 11 show that, unlike for intensity, the polarization asymmetry does not average out. Indeed, for the largest observer inclination (60°) shown in Figure 11, we see that there is a strong dependency of the degree of net polarization on the disk opening angle β. A similar but weaker effect is seen also for the observer inclinations i = 10° and i = 35°. Figure 11 also shows the behavior of the net polarization angle ψ. With all observer inclinations, a similar behavior of rotation of the polarization angle at high photon energies is seen. However, for these model parameters, the rotation mainly occurs at the high energy end of the spectrum, where the exponential cutoff makes the effect hard to observe in practice.

The changes in polarization with observation frequency described above, for β ∼ 0, are consistent with those of Schnittman & Krolik (2009), who studied an infinitely thin disk using a Monte Carlo approach. However, for a physically more realistic result, the accretion disk corona, as well as the radiation returning to and reflecting from the disk, needs to be taken into account, as in Schnittman & Krolik (2010). In addition, here we have shown that the geometry of the optically thick part of the disk cannot be ignored, which was an assumption in Schnittman & Krolik (2010). Combining the effects of the geometry with the effects of the corona and the returning radiation is straightforward using Arcmancer and will be investigated in a future work.

Another natural application of user-definable surfaces is the imaging of neutron stars. A solid surface is an excellent approximation for the radiating atmosphere of a neutron star, since the atmospheric thickness is on the order of ∼10 cm, whereas the radii of the neutron stars are in the range of ∼10 km (see e.g., Potekhin 2014 for a review). Thus, terminating geodesics on the top of the atmosphere and using a separate atmospheric model to provide the (angle-dependent) specific intensity and polarization as initial conditions are attractive possibilities.

The use of a numerical geodesic propagation code such as Arcmancer is further warranted due to the fact that a rotating neutron star is not exactly spherical but oblate, and the spacetime near the star cannot be exactly described by the Kerr metric (Stergioulas 2003; Bradley & Fodor 2009; Urbanec et al. 2013). Both complications are difficult to take into account when using fully analytic approaches, such as in, e.g., Pechenick et al. (1983), Strohmayer (1992), Miller & Lamb (1998), Poutanen & Gierliński (2003), and Lamb et al. (2009b), where the neutron star is modeled as a spherical surface in a Schwarzschild spacetime. The reason is twofold: the intersections of geodesics with the oblate surface are much more involved to compute (but not impossible; see Morsink et al. 2007; Lo et al. 2013; Miller & Lamb 2015; Stevens et al. 2016), and since Carter's constant (Carter 1968) of the Kerr solution is not available, the geodesics themselves cannot be analytically solved even in quadrature. Another benefit of using a fully covariant approach throughout is that the pitfalls of trying to combine special relativistic and general relativistic effects separately in an ad hoc way (as done in e.g., Lo et al. 2013) are avoided. For example, see Nättilä & Pihajoki (2018) and Lo et al. (2018) for a thorough discussion of an error in the calculation of the observed flux in the ad hoc approach that has gone undetected for years. Finally, incorporating polarization in an analytic geodesic propagator is only possible for Kerr (and Schwarzschild) spacetimes, but even then it is not trivial (see Viironen & Poutanen 2004; Dexter 2016). However, polarization data for this application are critical, since for small hot spots there is a severe degeneracy in the unpolarized pulse profile between the spot colatitude ${\theta }_{s}$ and the observer inclination i (Poutanen & Gierliński 2003).

In this section, we use Arcmancer to assess the effects of the oblateness of the neutron star surface and the deviation of the neutron star spacetime from the simple Schwarzschild spacetime on the radiative transfer calculation. For this purpose, we use the AlGendy–Morsink (AGM) form of the Butterworth–Ipser spacetime (AlGendy & Morsink 2014; see also Appendix C.3). The AGM spacetime describes the surroundings of a rotating neutron star, taking into account the oblate shape of the star. The spacetime is parametrized by the dimensionless rotational parameter $\bar{{\rm{\Omega }}}={\rm{\Omega }}{R}_{e}^{3/2}{M}^{-1/2}$ and the compactness parameter $x=M/{R}_{e}$, where Ω is the angular velocity of the rotation as seen by a distant observer, and M and R_e are the mass and the equatorial radius of the star, respectively. The oblate shape of the star is obtained from Equation (69).

As an example case, we studied a rotating neutron star with mass $M=1.6\,{M}_{\odot }$, equatorial radius ${R}_{e}=12\,\mathrm{km}$, and rotational frequency of $\nu =700\,\mathrm{Hz}$, with ${\rm{\Omega }}=2\pi \nu$. The high value of the spin was chosen to accentuate the effects of oblateness, yielding from Equation (70) a flattening of $f=1-{R}_{p}/{R}_{e}\sim 0.09$, where R_p and R_e are the polar and equatorial radii of the star, respectively. However, the high spin value is still within the observed range for neutron stars (Hessels et al. 2006). Similarly, the mass and the radius are well within the observed and inferred limits (Özel & Freire 2016; Steiner et al. 2016; Alsing et al. 2018). Using Arcmancer, we computed surface maps of the flux and polarization characteristics, again assuming that the emission originates from an electron scattering atmosphere, using Equations (21) and (22). This is a good approximation for neutron stars where the emission originates from thermonuclear outbursts on the surface (see e.g., Suleimanov et al. 2011 and the references therein). However, we note that the results can be extrapolated on a more qualitative level to shock-heated accretion-powered hot spots as well (see, e.g., Basko & Sunyaev 1976; Lyubarskii & Syunyaev 1982; Viironen & Poutanen 2004). Otherwise, the radiation transfer is computed as in Section 6.1.

The computations were repeated three times: for an oblate star using the AGM spacetime (hereafter, AGM+Obl), for an oblate star using the Schwarzschild spacetime (Sch+Obl), and for a spherical star in the Schwarzschild spacetime (Sch+Sph). The results are shown in Figures 12–15. Figure 12 shows the behavior of ${I}_{\nu }/{S}_{\nu }$, the specific intensity divided by the source function, and the polarization over the star surface, computed using the AGM spacetime at observer inclinations of $i=20^\circ$, 50°, and 90°. The combination of Doppler boosting and the strong angular dependence of the electron scattering atmosphere yields an intensity that varies significantly over the neutron star surface. The net polarization is high only near the edges, where the impact angle is large. Figure 12 also shows two possible paths of constant colatitude hot spots, assuming that the star is rotating around the vertical axis. From the figure, it is then easy to appreciate that a rotating hot spot should exhibit a large periodic variation in the observed polarization angle. This variation can be directly seen in Figure 15, which is consistent with the results in Viironen & Poutanen (2004).

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

Figure 13 shows the difference in normalized intensity ${I}_{\nu }/{S}_{\nu }$, degree of polarization, and polarization angle when the computation is performed using the AGM metric versus the Schwarzschild metric (i.e., AGM+Obl versus Sch+Obl). The effects of the rotation become significant only near the star, and consequently, the differences stay moderate for the most part, below ∼10% for the intensity and below ∼2% for the degree of polarization. There are areas of larger differences, but these are concentrated on the edges of the visible disk of the neutron star, and their total area is small. The differences in the polarization angle are larger, around ∼20° overall. There are very large differences near the point where the radiation was emitted toward the zenith in the frame of the neutron star surface, but this area corresponds to vanishing polarization, and as such, these differences are unobservable.

In contrast, Figure 14 displays the same differences but between computations performed using an oblate star versus a spherical star, both in Schwarzschild spacetimes (i.e., Sch+Obl versus Sch+Sph). The spherical star was given a radius equal to the equatorial radius of the oblate star. In this case, the differences in all quantities are much more pronounced. This is not a surprise, since a change in the shape of the star affects the redshift distribution on the surface due to variations in local surface gravity. These differences become even more evident when one looks at Figure 15, which shows two examples of light curves and the time-varying degree of polarization and polarization angle for a rotating hot spot. First, the pulse and polarization profiles closely match those obtained by Viironen & Poutanen (2004) for the Sch+Sph case and confirm that the observational degeneracy in unpolarized flux between observer inclination and spot colatitude is lifted by the polarization measurements. However, from the figure, it can be seen that the approximation of a spherical star produces results that differ significantly both in intensity and polarization properties from the result obtained using an oblate surface. In addition, there is a small but non-negligible difference between the results obtained using the AGM metric versus a plain Schwarzschild metric. Similar results for the unpolarized flux were obtained already in Psaltis & Özel (2014), although for an isotropically emitting atmosphere.

Based on our preliminary study, we can conclude that the error introduced when computing the polarization angle with the Schwarzschild spacetime approximation is largest when both the observer inclination and the spot colatitude are small. Likewise, the error in the degree of polarization is largest when the spot is near the equator, i.e., spot colatitude is close to ∼90°. We conclude that to obtain polarized pulse profiles that are accurate to below the ∼1% level, it is necessary that the rotation and the geometric shape of the star are both accurately modeled. In practice this means that the analytic results based on the Schwarzschild spacetime such as in, e.g., Weinberg et al. (2001), Viironen & Poutanen (2004), Lamb et al. (2009a), Lo et al. (2013), and Miller & Lamb (2015) should be used with caution. However, to actually reach the ∼1% level of accuracy, other systematic errors in, e.g., modeling the emission from the neutron star and its surrounding environment would also need to be resolved.

To further explore the possibility of using arbitrary metrics and multiple surfaces that may also move, we consider a toy model of an accreting black hole with a secondary black hole companion. To set up the problem, we use an approximative metric, constructed using the outgoing KS form of the Kerr metric, Equations (57) and (60). In the limit of zero spin, a = 0, the metric is

$$\begin{eqnarray}&&{g}_{{ab}}={\eta }_{{ab}}+F(M,r){l}_{a}({\boldsymbol{x}}){l}_{b}({\boldsymbol{x}}),\end{eqnarray} \tag{ 25 }$$

where M is the mass of the black hole, ${\boldsymbol{x}}=(x,y,z)$,

$$\begin{eqnarray}&&F(M,r)=-\displaystyle \frac{2M}{r},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{l}_{a}({\boldsymbol{x}})=\left(-1,\displaystyle \frac{x}{r},\displaystyle \frac{y}{r},\displaystyle \frac{z}{r}\right),\end{eqnarray} \tag{ 27 }$$

and ${r}^{2}={x}^{2}+{y}^{2}+{z}^{2}$. To this form we add a perturbation representing a distant second black hole moving at a slow coordinate velocity. Taking the spatial position of the second black hole to be a function ${{\boldsymbol{x}}}_{2}(t)$ of the coordinate time, we set

$$\begin{eqnarray}{g}_{{ab}} & = & {\eta }_{{ab}}+F({M}_{1},r){l}_{a}({\boldsymbol{x}}){l}_{b}({\boldsymbol{x}})\\ & & +\,F({M}_{2},{r}_{2}){l}_{a}({\boldsymbol{x}}-{{\boldsymbol{x}}}_{2}){l}_{b}({\boldsymbol{x}}-{{\boldsymbol{x}}}_{2}),\end{eqnarray} \tag{ 28 }$$

where M₁ and M₂ are the masses of the primary and secondary black holes, respectively, ${{\boldsymbol{x}}}_{2}$ is the spatial position of the secondary black hole, and ${r}_{2}^{2}={(x-{x}_{2})}^{2}$ $+\,{(y-{y}_{2})}^{2}\,+{(z-{z}_{2})}^{2}$.

The metric (28) is not a solution of the vacuum Einstein field equations, for which no exact dynamic binary black hole solution is known.⁶ For example, the metric (28) does not contain the gravitational wave component expected from the motion of multiple gravitating bodies. However, in the limit ${M}_{2}\sim 0$ and $d{{\boldsymbol{x}}}_{2}(t)/{dt}\sim 0$ for all t, the perturbation caused by the secondary is small and remains small, and the gravitational wave component is negligible, and in this sense the approximation is reasonable. For black hole binary systems with smaller separations and larger velocities, a discretized metric from a full GR simulation should be used with Arcmancer. More accurate analytical approximations, such as from Mundim et al. (2014), can also yield satisfactory accuracy (Sadiq et al. 2018), but the approximative analytical metrics are on the other hand algebraically complex.

We set ${M}_{1}=5\times {10}^{6}\,{M}_{\odot }$, representing an SMBH, and ${M}_{2}=0.05{M}_{1}$, which falls in the intermediate mass black hole (IMBH) range. Otherwise, we set up the system as in Section 6.1, by placing an image plane with physical dimensions ${[-50{M}_{1},50{M}_{1}]}^{2}$ at $R={10}^{7}{M}_{1}$ at an observer inclination of $i=60^\circ$. The secondary black hole is set on a rectilinear coordinate path ${{\boldsymbol{x}}}_{2}(t)={{\boldsymbol{x}}}_{\mathrm{2,0}}+t{{\boldsymbol{v}}}_{2}$, where

$$\begin{eqnarray}{{\boldsymbol{x}}}_{\mathrm{2,0}} & = & (L\sin i-\delta \cos i\cos {\phi }_{v},-\delta \sin {\phi }_{v},\\ & & \times \,L\cos i+\delta \sin i\cos {\phi }_{v}),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{2}={v}_{2}(-\cos i\sin {\phi }_{v},\cos {\phi }_{v},\sin i\sin {\phi }_{v}).\end{eqnarray} \tag{ 30 }$$

Here, δ is the apparent offset of the secondary's path, $L\,=\sqrt{{r}_{2,0}^{2}-{\delta }^{2}}$ is the orthogonal distance from the primary to the image plane, ${r}_{\mathrm{2,0}}$ is the minimum distance between the black holes, v₂ is the velocity of the secondary, and ${\phi }_{v}$ is the angle between the path of the secondary and the image plane x-axis. For this particular example, we set ${r}_{\mathrm{2,0}}={10}^{3}{M}_{1}$, $\delta \,=-4{M}_{1}$, ${\phi }_{v}=25^\circ$, and $v=3.162\times {10}^{-2}$. These initial conditions approximately correspond to an IMBH on a circular orbit around an SMBH, a situation that could possibly follow a merger of a more massive galaxy with a dwarf galaxy (Graham & Scott 2013). In order to have something to make a mock observation of, the primary black hole was given an infinitely thin Novikov–Thorne accretion disk, with $\alpha =0.1$ and an accretion rate $\dot{m}=0.01$ in units of the Eddington accretion rate. Geodesics were then propagated backwards in time from the image plane starting at 150 different values of the coordinate time, evenly distributed in $[-3000{M}_{1},3000{M}_{1}]$. For each set of geodesics, mock images, integrated fluxes, and polarization fraction and angle were computed.

The resulting light curves are shown in Figure 17, with selected resolved frames shown in Figure 16. The main effect of the passing secondary is a strong enhancement by a factor of ∼2 of the observed flux from the accretion disk of the primary, caused by gravitational lensing. The flux curve has a clearly non-sinusoidal shape, where, after the main peak, there is a pronounced shoulder. The double-peaked structure results from the lensing of the two main visible arcs of the primary accretion disk. The difference is that the major peak has a larger contribution from the Doppler-boosted side of the primary accretion disk. This asymmetry is also clearly visible in the curves for polarization fraction P and angle ψ; see Equations (23) and (24). The polarization fraction curve shows a clear two-peaked shape, with a sharp peak followed by a sharp trough. The polarization angle mirrors this behavior, with a maximum rotation of ∼7°.

Zoom In Zoom Out Reset image size

The light curves were also computed with a smaller value of observer inclination of i = 5°, also shown in Figure 17. The results show that the double-peaked structure of the light curve is more evident toward i = 0°, whereas the relative amount of polarized flux becomes significantly smaller. Both effects are to be expected considering the increased symmetry when $i\to {0}^{\circ }$. The changes in the degree of polarization and the polarization angle are more pronounced as well, but due to the negligible relative amount of polarized flux, these are unlikely to be detectable.

Zoom In Zoom Out Reset image size

Over longer timescales, the recurrent lensing by the secondary produces a periodic signal, which can be clearly observable over the baseline brightness of the primary accretion disk, as seen from Figure 17. However, the signal is strongly non-sinusoidal, which may reduce observability in periodicity searches based on periodogram techniques. On the other hand, if a series of accretion disk lensing events was observed, it should be possible to use lensing mock observation simulations to obtain independent constraints on the secondary black hole mass and the orbital parameters.

Finally, we note the interesting fact that the double-peaked light curve is reminiscent of the light curve of the periodic binary blazar OJ 287, which exhibits a long succession of strongly non-sinusoidal double-peaked outbursts every ∼12 years (Sillanpaa et al. 1988; Valtonen et al. 2008). Many different physical mechanisms for the outbursts have been proposed, such as tidally enhanced accretion rate (Sillanpaa et al. 1988), accretion disk impacts (Lehto & Valtonen 1996; Pihajoki 2016), and changes in the relativistic jet geometry (Katz 1997; Villata et al. 1998). Accretion disk lensing adds yet another possible outburst mechanism.

The basic building block of differential geometry is the manifold M, which can be intuitively understood as a space which locally "looks like" ${{\mathbb{R}}}^{n}$, the n-dimensional Euclidean space. More concretely, each manifold comes with an atlas of charts (coordinate systems) ${\phi }_{i}:M\supset {U}_{i}\to {{\mathbb{R}}}^{n}$, defined on open sets U_i of M. Using a chart ϕ, an abstract point $p\in M$ is transformed into its coordinate representation $({x}^{1}(p),\,\ldots ,\,{x}^{n}(p))\in {{\mathbb{R}}}^{n}$, where ${x}^{j}={\pi }^{j}\circ \phi$ is the projection to the jth coordinate. We say that the dimension of M is $\dim (M)=n$ and use the shorthand $\phi =({x}^{1},\,\ldots ,\,{x}^{n})$. A change in coordinates then corresponds to the transition map ${\phi }_{i}\circ {\phi }_{j}^{-1}:{{\mathbb{R}}}^{n}\to {{\mathbb{R}}}^{n}$, which changes a tuple of n coordinates to another tuple of n coordinates describing the same point. Arcmancer requires the manifold M to be differentiable, meaning that the maps ${\phi }_{i}\circ {\phi }_{j}^{-1}$ are differentiable. However, in the following we assume smooth (infinitely differentiable) manifolds.

The tangent space of M at point p, written T_pM, is the space of all vectors tangent to M at p. If $\dim (M)=n$, then T_pM is an n-dimensional real vector space. The dual space of T_pM, written ${T}_{p}^{* }M$ and called the cotangent space of M at p, is the space of all linear maps ${\omega }_{p}:{T}_{p}M\to {\mathbb{R}}$. These linear maps are called one-forms. For each chart $\phi =({x}^{1},\,\ldots ,\,{x}^{n})$, at point p we can define the vectors ${{\partial }_{i}| }_{p}:= {\left.\tfrac{\partial }{\partial {x}^{i}}\right|}_{p}$, called coordinate vectors. The vector ${{\partial }_{i}| }_{p}$ points toward the direction where the ith coordinate increases at p. The complete set of coordinate vectors, $\{{{\partial }_{1}| }_{p},\,\ldots ,\,{{\partial }_{n}| }_{p}\}$, forms a basis for T_pM, called the coordinate basis. Any vector $v\in {T}_{p}M$ can be written in terms of its components in this special basis as $v={v}^{i}{{\partial }_{i}| }_{p}$, where now ${v}^{i}\in {\mathbb{R}}$ are the components of v in the chart ϕ and the Einstein summation convention has been assumed. Likewise, there exists a set of one-forms, $\{{{{dx}}^{1}| }_{p},\,\ldots ,\,{{{dx}}^{n}| }_{p}\}$, called the coordinate one-forms, which forms a basis for ${T}_{p}^{* }M$ and makes it possible to write any one-form $\theta \in {T}_{p}^{* }M$ in terms of its components as $\theta ={\theta }_{i}{{dx}}^{i}$. The coordinate vectors and one-forms obey ${{{dx}}^{i}| }_{p}({{\partial }_{j}| }_{p})={\delta }_{j}^{i}$, where ${\delta }_{j}^{i}$ is the Kronecker delta symbol. Together, the coordinate vectors and one-forms are called the coordinate frame. The tangent and cotangent spaces can naturally have bases other than the coordinate bases as well (see Appendix A.6), but the coordinate basis is the default basis used in Arcmancer.

A rank (k, l) tensor can then be defined as a multilinear map $T:\mathop{\overbrace{{T}_{p}^{* }M\times \cdots \times {T}_{p}^{* }M}}\limits^{k\,\mathrm{times}}$ × $\mathop{\overbrace{{T}_{p}M\times \cdots \times {T}_{p}M}}\limits^{l\,\mathrm{times}}\to {\mathbb{R}}$. In general, the factors T_pM and ${T}_{p}^{* }M$ may be in any order, and the ordering is important, but for convenience, we will use the ordering given above. A rank $(0,0)$ corresponds to a scalar quantity, rank $(1,0)$ tensors are equivalent to vectors, and rank $(0,1)$ tensors to one-forms. The components of a tensor at p in the coordinate basis can be directly found from

$$\begin{eqnarray}&&{T}_{{i}_{k+1}\cdots +{i}_{k+l}}^{{i}_{1}\cdots {i}_{k}}=T({{{dx}}^{{i}_{k+1}}| }_{p},\,\ldots \,,{{{dx}}^{{i}_{k+l}}| }_{p},{{\partial }_{{i}_{1}}| }_{p},\,\ldots \,,{{\partial }_{{i}_{k}}| }_{p}).\end{eqnarray} \tag{ 31 }$$

Using the components, a tensor can be locally defined as an expansion

$$\begin{eqnarray}T & = & {T}_{{i}_{k+1}\cdots +{i}_{k+l}}^{{i}_{1}\cdots {i}_{k}}{{\partial }_{{i}_{1}}| }_{p}\otimes \cdots \otimes {{\partial }_{{i}_{k}}| }_{p}\\ & & \otimes \,{{{dx}}^{{i}_{k+1}}| }_{p}\otimes \cdots \otimes {{{dx}}^{{i}_{k+l}}| }_{p}.\end{eqnarray} \tag{ 32 }$$

If instead of the chart $\phi =({x}^{1},\,\ldots ,\,{x}^{n})$ we wish to use another (overlapping) chart $\psi =({y}^{1},\,\ldots ,\,{y}^{n})$ to represent the tensor T, the components ${T}_{{i}_{k+1}\cdots +{i}_{k+l}}^{{i}_{1}\cdots {i}_{k}}\in {\mathbb{R}}$ of a tensor must be transformed accordingly. The new components turn out to be

$$\begin{eqnarray}&&{T}_{{i}_{k+1}\cdots +{i}_{k+l}}^{{i}_{1}\cdots {i}_{k}}={T}_{{j}_{k+1}\cdots +{j}_{k+l}}^{{j}_{1}\cdots {j}_{k}}{J}_{{j}_{1}}^{{i}_{1}}\cdots {J}_{{j}_{k}}^{{i}_{k}}{({J}^{-1})}_{{i}_{k+1}}^{{j}_{k+1}}\cdots {({J}^{-1})}_{{i}_{k+l}}^{{j}_{k+l}},\end{eqnarray} \tag{ 33 }$$

where ${J}_{j}^{i}=\displaystyle {\rm{\partial }}({y}^{i}\circ {\phi }^{-1})/{\rm{\partial }}{x}^{j}$ are the components of the Jacobian J of the function $\psi \circ {\phi }^{-1}$ at point p, and similarly for the inverse of the Jacobian ${J}^{-1}$.

The definitions above generalize to vector, one-form, and tensor fields, which can be understood as functions which, for each point p of M, pick a specific vector, one-form, or tensor, respectively. In physics-oriented GR literature, all tensorial quantities are usually tensor fields. In addition, a convention called the abstract index notation (Penrose & Rindler 1987) is often used. In this convention, for example, a rank $(2,3)$ tensor field T can be written as $T{{}_{{cd}}^{{ab}}}^{e}$. The number and ordering of the upper and lower indexes is taken only to signify the number and ordering of the factors T_pM and ${T}_{p}^{* }M$ in the definition, Equation (32), of the tensor. This approach makes it possible to write all coordinate-invariant tensor operations tersely, without specifying any underlying basis. In this paper, we use the abstract index notation wherever possible.

A manifold may have a special rank $(0,2)$ tensor field called the metric, usually written g_ab. The metric defines the inner product of vectors $\langle {v}^{a}\,{w}^{a}\rangle ={g}_{{ab}}{v}^{a}{w}^{b}$ and consequently a norm $\parallel {v}^{a}\parallel =\sqrt{| {g}_{{ab}}{v}^{a}{v}^{b}| }$ on each T_pM. Intuitively, the metric defines the distance between nearby points x and $x+{\rm{\Delta }}x$ as the norm of the tangent vector approximated by ${\rm{\Delta }}x$. In physics, the components of the metric are often written in the form of a line element ${{ds}}^{2}$, essentially an expansion in terms of coordinate basis tensors, as ${{ds}}^{2}={g}_{{ij}}\,{{dx}}^{i}\otimes {{dx}}^{j}$. The signature of the metric is the pair (p, q) of the number of positive and negative eigenvalues of the matrix of components of g_ab (in any basis), respectively, assuming $p+q=n$. If q = 0, the metric and the manifold are said to be Riemannian. For signatures of $(n-1,1)$ or $(1,n-1)$, the metric and the manifold are said to be Lorentzian. In particular, GR is defined in terms of a four-dimensional Lorentzian manifold. In this paper, Lorentzian metrics are always assumed to be of type $(1,n-1)$, and in particular for GR, this implies the $(+---)$ metric convention.

A choice of metric also defines a unique metric-compatible zero-torsion connection ∇, often called the Levi–Civita connection. A connection in general can be roughly said to characterize which vectors of the tangent spaces of two nearby points are to be considered equal, or alternatively, to yield a vector field ${{\rm{\nabla }}}_{X}Y$ describing the rate of change of the vector field Y in the direction of the vector field X, called the (natural) covariant derivative of Y with respect to X. Written in terms of the coordinate vector fields, the connection is ${{\rm{\nabla }}}_{{\partial }_{i}}({\partial }_{j})={{\rm{\Gamma }}}_{{ij}}^{k}{\partial }_{k}$, where ${{\rm{\Gamma }}}_{{ij}}^{k}$ are the Christoffel symbols (of the second kind). This work only uses the Levi–Civita connection, the zero-torsion property of which can be written as ${{\rm{\Gamma }}}_{{ij}}^{k}={{\rm{\Gamma }}}_{{ji}}^{k}$, and the metric compatibility as ${{\rm{\nabla }}}_{c}{g}_{{ab}}=0$. In terms of the metric, the Christoffel symbols read

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}^{k}=\displaystyle \frac{1}{2}{g}^{{km}}({\partial }_{i}{g}_{{jm}}+{\partial }_{j}{g}_{{im}}-{\partial }_{m}{g}_{{ij}}).\end{eqnarray} \tag{ 34 }$$

The covariant derivative can be generalized for any tensor field T. Using the Christoffel symbols, the covariant derivative of T with respect to X can be written in component form as

$$\begin{eqnarray}{{\rm{\nabla }}}_{X}T & = & {X}^{i}{{\rm{\nabla }}}_{{\partial }_{i}}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}\\ & = & {X}^{i}({\partial }_{i}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}+{T}_{{b}_{1}\cdots {b}_{l}}^{{c}_{1}\cdots {a}_{k}}{{\rm{\Gamma }}}_{{c}_{1}i}^{{a}_{1}}\\ & & +\,\cdots +{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {c}_{k}}{{\rm{\Gamma }}}_{{c}_{k}i}^{{a}_{k}}-{T}_{{c}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}{{\rm{\Gamma }}}_{{b}_{1}i}^{{c}_{1}}\\ & & -\,\cdots -{T}_{{b}_{1}\cdots {c}_{l}}^{{a}_{1}\cdots {a}_{k}}{{\rm{\Gamma }}}_{{b}_{l}i}^{{c}_{l}}).\end{eqnarray} \tag{ 35 }$$

In the abstract index notation, the covariant derivative is written simply as ${{\rm{\nabla }}}_{X}T={X}^{c}{{\rm{\nabla }}}_{c}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}$. The term covariant derivative of T, without additional qualifiers, is often used to refer only to the vector field independent part ${{\rm{\nabla }}}_{c}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}$. If ${{\rm{\nabla }}}_{X}T=0$, we say that T is parallel (transported) along the vector field X. This concept can be extended to the case where X is a tangent vector field of a curve γ; see below.

For semi-Riemannian spaces, the metric also divides vectors at a point into three categories. A vector $v\in {T}_{p}M$ is said to be timelike if the inner product ${g}_{{ab}}{v}^{a}{v}^{b}$ is positive, null if it is zero, and spacelike otherwise. The same classification can be extended to vector fields V if the sign of the inner product is the same everywhere on M.

For ray-tracing mock observations, the notion of curves and especially geodesics on a manifold is essential. A curve γ can be thought of as a map $\gamma :{\mathbb{R}}\supset I\to M$. The curve defines a vector field u^a through $u(\lambda )=d\gamma (\lambda )/d\lambda$ along the curve. This vector field can be understood as the velocity vector field of the curve. The curve itself is a solution of

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}^{a}(\lambda )}{d\lambda }={{\rm{\nabla }}}_{{u}^{b}}{u}^{a}={u}^{b}{{\rm{\nabla }}}_{b}{u}^{a}.\end{eqnarray} \tag{ 36 }$$

If in particular ${u}^{b}{{\rm{\nabla }}}_{b}{u}^{a}=0$, the tangent vector of the curve is parallel transported along the curve, and we say that the curve is a geodesic. Equation (36) can also be written as

$$\begin{eqnarray}&&{u}^{b}{{\rm{\nabla }}}_{b}{u}^{a}={f}^{a},\end{eqnarray} \tag{ 37 }$$

where the vector field f^a is analogous to a force that causes the curve to deviate from a straight path along the manifold. For geodesics, f^a = 0, corresponding to the notion of geodesics as the straightest possible paths. Written in terms of components in a specific chart, Equation (37) reads

$$\begin{eqnarray}&&\displaystyle \frac{{{du}}^{a}(\lambda )}{d\lambda }=-{{\rm{\Gamma }}}_{{bc}}^{a}{u}^{b}{u}^{c}+{f}^{a}.\end{eqnarray} \tag{ 38 }$$

In this form, the equation can be solved numerically, as long as the solution stays within the domain of the chart.

Generic tensorial quantities can also be parallel transported along the curve. For a rank k + l tensor field ${T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}$, the equation of parallel transport is

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{{u}^{c}}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}={u}^{c}{{\rm{\nabla }}}_{c}{T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}=0,\end{eqnarray} \tag{ 39 }$$

which in component form can be found from Equation (35) with the substitution $X={u}^{a}$. Equation (39) can likewise be directly solved numerically.

Finally, curves for which the tangent vector field u^a is always timelike, spacelike, or null, are called timelike, spacelike, or null, respectively. This characterization is important for geodesics, for which the tangent vector field can never change their timelike, spacelike or null character. As such, geodesics always fall in one of these categories.

Manifolds may contain many kinds of submanifolds. A particularly useful class of submanifolds are those defined by level sets of functions, called level hypersurfaces. If $S:M\to {\mathbb{R}}$ is a smooth map, then the sets ${S}_{c}:= {S}^{-1}(c)=\{p\in M| S(p)=c\}$ are submanifolds of M for each c in the codomain of S if ${{dS}| }_{p}\ne 0$ for all $p\in {S}_{c}$. We can subsume the constant c in the definition of the function S and take c = 0, which is always assumed in the Arcmancer code. Furthermore, with a slight abuse of notation, we use the function S defining the hypersurface to refer to the hypersurface itself.

A level hypersurface divides the manifold M into two disconnected subsets corresponding to regions where the value of S is either negative or positive. In some cases, these may be conveniently taken to be the "outside" and the "inside" of a region bounded by the level hypersurface.

Tensor fields can also be expressed in terms of bases other than the coordinate basis. A local frame defined in an open set $U\supset M$ consists of n vector fields $\{{E}_{1},\,\ldots ,\,{E}_{n}\}$ that form a basis for T_pM at each $p\in U$. For each local frame, there is a corresponding local coframe of one-form fields $\{{\omega }^{1},\,\ldots ,\,{\omega }^{n}\}$ for which ${\omega }^{i}({E}_{j})={\delta }_{i}^{j}$. The components of a rank (k, l) tensor ${T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}$ in terms of these bases can then be defined through (see Equation (31))

$$\begin{eqnarray}&&{T}_{{j}_{1}\cdots {j}_{l}}^{{i}_{1}\cdots {i}_{k}}=T({\omega }^{{i}_{1}},\,\ldots ,\,{\omega }^{{i}_{k}},{E}_{{j}_{1}},\,\ldots ,\,{E}_{{j}_{l}}).\end{eqnarray} \tag{ 40 }$$

The transformation from the components of a tensor in a coordinate frame to components in terms of a local frame can be found by first writing the local (co)frame in terms of the coordinate frame as ${E}_{i}={\sum }_{\alpha }{E}_{i}^{\alpha }{\partial }_{\alpha }$ and ${\omega }^{j}={\sum }_{\beta }{\omega }_{\beta }^{j}{{dx}}^{\beta }$. Using the multilinearity of an arbitrary (k, l) tensor ${T}_{{b}_{1}\cdots {b}_{l}}^{{a}_{1}\cdots {a}_{k}}$ in Equation (40) then yields

$$\begin{eqnarray}{T}_{{j}_{1}\cdots {j}_{l}}^{{i}_{1}\cdots {i}_{k}} & = & \displaystyle \sum _{{\beta }_{1}}\cdots \displaystyle \sum _{{\alpha }_{l}}{\omega }_{{\beta }_{1}}^{{i}_{1}}\cdots {\omega }_{{\beta }_{k}}^{{i}_{k}}{E}_{{j}_{1}}^{{\alpha }_{1}}\cdots {E}_{{j}_{l}}^{{\alpha }_{l}}\\ & & \times \,T({{dx}}^{{\beta }_{1}},\,\ldots ,\,{{dx}}^{{\beta }_{k}},{\partial }_{{\alpha }_{1}},\,\ldots ,\,{\partial }_{{\alpha }_{l}}),\end{eqnarray} \tag{ 41 }$$

where ${T}_{{\alpha }_{1}\cdots {\alpha }_{l}}^{{\beta }^{1}\cdots {\beta }^{k}}=T({{dx}}^{{\beta }_{1}},\,\ldots ,\,{{dx}}^{{\beta }_{k}},{\partial }_{{\alpha }_{1}},\,\ldots ,\,{\partial }_{{\alpha }_{l}})$ are the components of the tensor in the coordinate frame. From a computational perspective, it is useful to note that the set of numbers ${E}_{i}^{\alpha }$ can be interpreted as a matrix, for which the corresponding set ${\omega }_{\beta }^{j}$ forms an inverse matrix. Equation (41) is seen to resemble Equation (33), with these matrices taking the place of the Jacobians.

Orthonormal local frames define a set of convenient projections of T_pM onto orthogonal complements at each point p. A collection of k non-null basis vectors ${ \mathcal K }\,=\{{E}_{{i}_{j}}\in {T}_{p}M| j=1,\,\ldots ,\,k\}$ defines a projection onto the subspace of T_pM spanned by ${ \mathcal K }$ through

$$\begin{eqnarray}&&{P}_{\parallel }{({E}_{{i}_{1}},\,\ldots ,\,{E}_{{i}_{k}})}_{b}^{a}=\displaystyle \sum _{j=1}^{k}\displaystyle \frac{{E}_{{i}_{j}}^{a}{({E}_{{i}_{j}})}_{b}}{{E}_{{i}_{j}}^{c}{({E}_{{i}_{j}})}_{c}}.\end{eqnarray} \tag{ 42 }$$

Similarly, a projection onto the orthogonal complement of this subspace is defined via

$$\begin{eqnarray}&&{P}_{\perp }{({E}_{{i}_{1}},\,\ldots ,\,{E}_{{i}_{k}})}_{b}^{a}={\delta }_{b}^{a}-\displaystyle \sum _{j=1}^{k}\displaystyle \frac{{E}_{{i}_{j}}^{a}{({E}_{{i}_{j}})}_{b}}{{E}_{{i}_{j}}^{c}{({E}_{{i}_{j}})}_{c}}.\end{eqnarray} \tag{ 43 }$$

For Lorentzian spaces, direct projections with respect to a null vector k^a do not work. However, given a unit timelike vector u^a, a projection onto the space orthogonal to both u^a and k^a can be given through

$$\begin{eqnarray}{P}_{\perp }{(u,k)}_{b}^{a} & = & {\delta }_{b}^{a}-\displaystyle \frac{{u}^{a}{u}_{b}}{{u}^{c}{u}_{c}}-\displaystyle \frac{{s}^{a}{s}_{b}}{{s}^{c}{s}_{c}}\\ & = & {\delta }_{b}^{a}-\displaystyle \frac{{u}^{a}{k}_{b}}{K}+\displaystyle \frac{{k}^{a}{u}_{b}}{K}+\displaystyle \frac{{k}^{a}{k}_{b}}{{K}^{2}},\end{eqnarray} \tag{ 44 }$$

where ${s}^{a}={P}_{\perp }(u{)}_{b}^{a}{k}^{b}$ and $K={u}_{a}{k}^{a}$. In four dimensions, the operator ${P}_{\perp }(u,k{)}_{b}^{a}$ is a projection onto a two-dimensional surface on which the observed components of polarization are defined. In this context, ${P}_{\perp }(u,k{)}_{b}^{a}$ is known as the screen projection operator (Gammie & Leung 2012).

A.6.1. Lorentz Frames and Observers

In the case of a four-dimensional Lorentzian manifold, an orthonormal local frame is often called a (local) Lorentz frame or a tetrad. The tetrad can in general consist of any permissible combination of null, spacelike, and timelike vector fields. However, in this work, we use the term Lorentz frame to mean a mutually orthonormal combination of one timelike vector, E_t, and three spacelike vectors, E_x, E_y, and E_z, defined at a point. Since parallel transport preserves inner products, a parallel transported Lorentz frame is still a valid Lorentz frame.

Observers in GR are characterized by a timelike curve γ, a world line. An observer's rest frame can then be defined as a Lorentz frame for which the timelike basis vector is given by the observer's four-velocity, or ${E}_{t}=d\gamma (\lambda )/d\lambda$, and the spatial triad $\{{E}_{x},{E}_{y},{E}_{z}\}$ defines the observer's choice of spatial coordinate system. The values of physical quantities as measured by the observer are obtained by expressing them in the observer's rest-frame basis.

For example, an important feature of GR is that angles between vectors at a point is observer dependent. The angle θ between vectors x^a and y^a as measured by an observer with a four-velocity u^a is then

$$\begin{eqnarray}&&\cos \theta (u;x,y)=\displaystyle \frac{-{x}_{\perp }^{a}{y}_{\perp a}}{\sqrt{| ({x}_{\perp }^{b}{x}_{\perp b})({y}_{\perp }^{c}{y}_{\perp c})| }},\end{eqnarray} \tag{ 45 }$$

where ${x}_{\perp }^{a}={P}_{\perp }(u{)}_{b}^{a}{x}^{b}$ and ${y}_{\perp }^{a}={P}_{\perp }(u{)}_{b}^{a}{y}^{b}$.

Briefly, a propagating monochromatic radiation front can be modeled as a congruence of curves, each perpendicular to the surface of constant phase. This is possible in the limit where the wavelength of the radiation is much smaller than the scale of variations in the radiation front (curvature, amplitude, polarization) and much smaller than the local "radius of curvature of the space" $\propto {| R{}_{{bcd}}^{a}| }^{-1/2}$, where $R{}_{{bcd}}^{a}$ is the Riemann curvature tensor.

Furthermore, if the contribution of matter is insignificant, i.e., ${T}^{{ab}}\sim 0$, the result is that the propagation of radiation can be modeled by solving an equation of radiative transfer along null geodesics. If the integrated contribution of the intervening matter is non-negligible, the radiation front normals will not be geodesics in general, and there will be a matter-dependent forcing term f^a in Equation (37) (see, e.g., Broderick & Blandford 2003, 2004).

The classical equation of radiative transfer in Cartesian coordinates is (Mihalas & Mihalas 1984)

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{c\,\partial t}+\hat{{\boldsymbol{n}}}\cdot {\rm{\nabla }}\right){I}_{\nu }({\boldsymbol{x}},t;\hat{{\boldsymbol{n}}})\\ &&\quad =\,{j}_{\nu }({\boldsymbol{x}},t;\hat{{\boldsymbol{n}}})-{\alpha }_{\nu }({\boldsymbol{x}},t;\hat{{\boldsymbol{n}}}){I}_{\nu }({\boldsymbol{x}},t;\hat{{\boldsymbol{n}}}),\end{eqnarray} \tag{ 46 }$$

where j_ν is the total emission coefficient (emissivity) and ${\alpha }_{\nu }$ the total absorption coefficient. Equation (46) describes the change in specific intensity at a frequency ν in the direction $\hat{{\boldsymbol{n}}}$ at a point ${\boldsymbol{x}}$ and time t. Often in astrophysical cases, the one-dimensional time-independent case suffices, in which case Equation (46) reduces to a more commonly seen form,

$$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{\nu }(s)}{{ds}}={j}_{\nu }(s)-{\alpha }_{\nu }(s){I}_{\nu }(s),\end{eqnarray} \tag{ 47 }$$

where s is the distance along a radiation front normal. Equations (46) and (47) do not take into account interference, quantum effects, or, most importantly, polarization.

Polarization can be included by making the specific intensity I_ν vector-valued by introducing the Stokes intensities ${{\boldsymbol{I}}}_{\nu }\,=({I}_{\nu },{Q}_{\nu },{U}_{\nu },{V}_{\nu })$. Similarly, j is replaced by the vector

$$\begin{eqnarray}&&{{\boldsymbol{J}}}_{\nu }=({j}_{I,\nu },{j}_{Q,\nu },{j}_{U,\nu },{j}_{V,\nu })\end{eqnarray} \tag{ 48 }$$

of Stokes emissivities and α by the response (or Müller) matrix

$$\begin{eqnarray}{{\boldsymbol{M}}}_{\nu }=\left(\begin{array}{cccc}{\alpha }_{I,\nu } & {\alpha }_{Q,\nu } & {\alpha }_{U,\nu } & {\alpha }_{V,\nu }\\ {\alpha }_{Q,\nu } & {\alpha }_{I,\nu } & {r}_{V,\nu } & -{r}_{U,\nu }\\ {\alpha }_{U,\nu } & -{r}_{V,\nu } & {\alpha }_{I,\nu } & {r}_{Q,\nu }\\ {\alpha }_{V,\nu } & {r}_{U,\nu } & -{r}_{Q,\nu } & {\alpha }_{I,\nu }\end{array}\right),\end{eqnarray} \tag{ 49 }$$

where the α-coefficients represent absorption effects, and the r-coefficients relate to Faraday conversion and rotation. It should be noted that different conventions for ${\boldsymbol{M}}$ exist, varying by the sign of r_U. The one-dimensional polarized equation of radiative transfer is then

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{I}}}_{\nu }(s)}{{ds}}={{\boldsymbol{J}}}_{\nu }-{{\boldsymbol{M}}}_{\nu }(s){{\boldsymbol{I}}}_{\nu }(s).\end{eqnarray} \tag{ 50 }$$

The general relativistic generalization of Equation (50) is (e.g., Gammie & Leung 2012)

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}^{{ab}}}{d\lambda }={J}^{{ab}}+{H}^{{abcd}}{N}_{{cd}},\end{eqnarray} \tag{ 51 }$$

where λ is the affine parameter along the curve representing the propagating radiation front, N^ab is the (complex-valued) polarization tensor, J^ab is the emissivity tensor, and H^abcd is the response tensor.

Directly solving Equation (51) requires integrating the 16 real independent components of N^ab. This number can be reduced to four by parallel transporting a polarization frame along the geodesic. The frame consists of two orthogonal spacelike vectors that are also orthogonal to the geodesic and the observer four-velocity. Expressing all quantities in this frame using the screen projection operator (44), Equation (51) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{ \mathcal I }}}_{\nu }}{d\lambda }={{\boldsymbol{ \mathcal J }}}_{\nu }-{{\boldsymbol{ \mathcal M }}}_{\nu }{{\boldsymbol{ \mathcal I }}}_{\nu },\end{eqnarray} \tag{ 52 }$$

where ${{\boldsymbol{ \mathcal I }}}_{\nu }={\nu }^{-3}{{\boldsymbol{I}}}_{\nu }$, ${\boldsymbol{ \mathcal J }}=C{\nu }^{-2}{{\boldsymbol{J}}}_{\nu }$, ${\boldsymbol{ \mathcal M }}=C{\nu }^{-1}{\boldsymbol{M}}$, and C is a constant related to the parametrization of the curve. An observer with a four-velocity u^a at one end of the curve, where the tangent is k^a, has $C={u}_{a}{k}^{a}/{\nu }_{0}$ for an observed frequency ${\nu }_{0}$.

C.1.1. Two-sphere

Spherical coordinates $(\theta ,\phi )$—Arcmancer uses two copies of the spherical chart $(\theta ,\phi )$ to cover the entire two-sphere. The metric, in either copy, is given by the usual

$$\begin{eqnarray}&&{{ds}}^{2}=d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}.\end{eqnarray} \tag{ 53 }$$

C.2.1. Minkowski Spacetime

Cartesian coordinates $(t,x,y,z)$— The Cartesian Minkowski coordinates are often in the literature denoted by ${\eta }_{{ab}}$, and we use the same convention. The line element is diagonal, given by

$$\begin{eqnarray}&&{{ds}}^{2}={{dt}}^{2}-{{dx}}^{2}-{{dy}}^{2}-{{dz}}^{2}.\end{eqnarray} \tag{ 54 }$$

Spherical coordinates $(t,r,\theta ,\phi )$—

$$\begin{eqnarray}&&{{ds}}^{2}={{dt}}^{2}-{{dr}}^{2}-{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta \,d{\phi }^{2})\end{eqnarray} \tag{ 55 }$$

C.2.2. Kerr Spacetime

The Kerr spacetime (Kerr 1963) is an important solution to the Einstein field equations, representing a rotating black hole exactly and other rotating fluid bodies of finite size asymptotically. The Kerr spacetime is parametrized by the mass M and the angular momentum J. Usually, J is given through the normalized spin parameter $a=J/M$, in which case $a\in [0,M]$, or through the dimensionless spin parameter χ, so that $\chi \in [0,1]$. The solution reduces to the Schwarzschild spacetime when χ = 0 and further to the flat Minkowski spacetime when M = 0.

Boyer–Lindquist coordinates $(t,r,\theta ,\phi )$—Perhaps the clearest representation of the Kerr metric in BL coordinates is through the matrix form

$$\begin{eqnarray}({g}_{ab})=\left(\begin{array}{cccc}1-\displaystyle \frac{2Mr}{{\rho }^{2}} & 0 & 0 & \displaystyle \frac{2Mra{\sin }^{2}\theta }{{\rho }^{2}}\\ 0 & -\displaystyle \frac{{\rho }^{2}}{{\rm{\Delta }}} & 0 & 0\\ 0 & 0 & -{\rho }^{2} & 0\\ \displaystyle \frac{2Mra{\sin }^{2}\theta }{{\rho }^{2}} & 0 & 0 & -{\sin }^{2}\theta \left({r}^{2}+{a}^{2}+\displaystyle \frac{2{Mra}^{2}{\sin }^{2}\theta }{{\rho }^{2}}\right)\,\end{array}\right),\end{eqnarray} \tag{ 56 }$$

where ${\rho }^{2}={r}^{2}+{a}^{2}{\cos }^{2}\theta$ and ${\rm{\Delta }}={r}^{2}-2{Mr}+{a}^{2}$. The BL form of the metric is singular when ${\rm{\Delta }}=0$ or ${\rho }^{2}=0$. The former condition corresponds to ${r}_{\pm }={m}^{2}\pm \sqrt{{m}^{2}-{a}^{2}}$, which gives the locations of the inner and outer event horizons, where the curvature is not singular. The condition ${\rho }^{2}=0$ implies r = 0 and $\cos \theta =0$, and corresponds to the curvature singularity.

Cartesian KS coordinates $(t,x,y,z)$—The Cartesian KS coordinates come in two flavors, ingoing and outgoing, adapted to null geodesics that move radially inwards or outwards, respectively. As such, they are a generalization of the ingoing and outgoing Eddington–Finkelstein coordinates for a Schwarzschild black hole. Of the two variants, only the ingoing form is typically seen in the literature. The metric in the ingoing KS coordinates is most easily written as a sum,

$$\begin{eqnarray}&&{g}_{{ab}}={\eta }_{{ab}}-F(x,y,z;a){l}_{a}{l}_{b},\end{eqnarray} \tag{ 57 }$$

where η is the Cartesian Minkowski metric,

$$\begin{eqnarray}&&F(x,y,z;a)=\displaystyle \frac{2{{Mr}}^{3}}{{r}^{4}+{a}^{2}{z}^{2}},\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&({l}_{a})=\left(1,\displaystyle \frac{{rx}+{ay}}{{a}^{2}+{r}^{2}},\displaystyle \frac{{ry}-{ax}}{{a}^{2}+{r}^{2}},\displaystyle \frac{z}{r}\right),\end{eqnarray} \tag{ 59 }$$

and r is defined implicitly through ${x}^{2}+{y}^{2}+{z}^{2}\,={r}^{2}+{a}^{2}(1-{z}^{2}/{r}^{2})$. The vector l_a is null with respect to both ${\eta }_{{ab}}$ and g_ab. The metric in the outgoing coordinates is similarly given as a sum as in Equation (57), but with

$$\begin{eqnarray}&&({l}_{a})=\left(-1,\displaystyle \frac{{rx}-{ay}}{{a}^{2}+{r}^{2}},\displaystyle \frac{{ry}+{ax}}{{a}^{2}+{r}^{2}},\displaystyle \frac{z}{r}\right).\end{eqnarray} \tag{ 60 }$$

When numerically computing geodesics near a Kerr black hole, it is crucial that the right form of the KS coordinates is chosen. If the geodesic approaches the hole, the ingoing chart should be used, and for a geodesic going away from the hole, the outgoing chart should be used. Failure to do so has significant computational penalties, as can be seen in Section 5.1.2. Note that this implies that for a geodesic coming in toward a Kerr black hole, passing by it, and then leaving, both ingoing and outgoing charts should be used with a change of coordinates near the closest approach.

AGM coordinates $(t,\bar{r},\theta ,\phi )$—The AGM spacetime and coordinates (AlGendy & Morsink 2014) actually refer to a specific choice of a Butterworth–Ipser (BI) spacetime, with the accompanying coordinate chart (Butterworth & Ipser 1976). The BI spacetime is a general representation of the spacetime outside an axisymmetric rotating fluid body. The AGM spacetime is a special case representing the spacetime around a physically realistic oblate rotating neutron star. The AGM representation is accurate up to second order in the dimensionless rotation parameter $\bar{{\rm{\Omega }}}={\rm{\Omega }}{R}_{e}^{3/2}{M}^{-1/2}$, where Ω is the rotational angular velocity of the star as seen by a distant observer, R_e is the equatorial radius of the star, and M is the mass of the star. In BI coordinates, the generic axisymmetric metric reads

$$\begin{eqnarray}({g}_{{ab}})=\left(\begin{array}{cccc}{e}^{2\nu }-{B}^{2}{\bar{r}}^{2}{\omega }^{2}{\sin }^{2}\theta {e}^{-2\nu } & 0 & 0 & {B}^{2}{\bar{r}}^{2}\omega {\sin }^{2}\theta {e}^{-2\nu }\\ 0 & -{e}^{2(\zeta -\nu )} & 0 & 0\\ 0 & 0 & -{e}^{2(\zeta -\nu )}{\bar{r}}^{2} & 0\\ {B}^{2}{\bar{r}}^{2}\omega {\sin }^{2}\theta {e}^{-2\nu } & 0 & 0 & -{B}^{2}{\bar{r}}^{2}{\sin }^{2}\theta {e}^{-2\nu }\end{array}\right),\end{eqnarray} \tag{ 61 }$$

where ν, B, ω, and ζ are the so-called metric functions or potentials. The AGM metric is defined by using

$$\begin{eqnarray}&&\nu =\mathrm{ln}\displaystyle \frac{1-M/(2\bar{r})}{1+M/(2\bar{r})}+\left(\displaystyle \frac{\beta }{3}-{{qP}}_{2}(\cos \theta )\right){\left(\displaystyle \frac{M}{\bar{r}}\right)}^{3},\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&B=1-{\left(\displaystyle \frac{M}{2\bar{r}}\right)}^{2}+\beta {\left(\displaystyle \frac{M}{\bar{r}}\right)}^{2},\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&\omega =\displaystyle \frac{2j}{\bar{r}}{\left(\displaystyle \frac{M}{\bar{r}}\right)}^{2}\left(1-\displaystyle \frac{3M}{\bar{r}}\right),\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&\zeta =\mathrm{ln}\left[1-{\left(\displaystyle \frac{M}{2\bar{r}}\right)}^{2}\right]+\beta \displaystyle \frac{4{P}_{2}(\cos \theta )-1}{3},\end{eqnarray} \tag{ 65 }$$

where P₂ is the second-order Legendre polynomial. The constants q and β are the dimensionless moments of energy density and pressure, respectively, and $j=J/{M}^{2}$ is the dimensionless angular momentum. AlGendy & Morsink (2014) found these constants to be well described across various neutron star parameters and equations of state by the approximate relations

$$\begin{eqnarray}&&j=(1.136-2.53x+5.6{x}^{2})\bar{{\rm{\Omega }}},\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&q=-0.11{x}^{-2}{\bar{{\rm{\Omega }}}}^{2},\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&\beta =0.4454x{\bar{{\rm{\Omega }}}}^{2},\end{eqnarray} \tag{ 68 }$$

where $x=M/{R}_{e}$ is called the compactness (parameter). Along with the metric potentials, AlGendy & Morsink (2014) also derived an equation for the shape of the surface of a rotating neutron star, used in Section 6.2,

$$\begin{eqnarray}&&R(\theta )={R}_{e}[1+(-0.788+1.030x){\bar{{\rm{\Omega }}}}^{2}{\cos }^{2}\theta ].\end{eqnarray} \tag{ 69 }$$

From this, the flattening, also called the oblateness, of the star is given by

$$\begin{eqnarray}&&f=\displaystyle \frac{R(\pi /2)-R(0)}{R(\pi /2)}=(0.788-1.030x){\bar{{\rm{\Omega }}}}^{2}.\end{eqnarray} \tag{ 70 }$$

It should be noted that the quantities R_e and $R(\theta )$ are not defined in terms of $\bar{r}$, but instead in terms of a radial coordinate $r={{Be}}^{-\nu }\bar{r}$.

The Hartle–Thorne spacetime (Hartle & Thorne 1968) describes the spacetime around a rotating oblate star. The original derivation was based on a perturbation of the non-rotating Schwarzschild spacetime. In Arcmancer, we use instead a version based on a perturbation of the rotating Kerr spacetime, given in Glampedakis & Babak (2006).

Glampedakis–Babak (GB) coordinates $(t,r,\theta ,\phi )$—The GB coordinates are based on the BL coordinates of the Kerr spacetime. In these coordinates, in a parametrization used by Bauböck et al. (2012), the metric is given by

$$\begin{eqnarray}&&{g}_{{ab}}={g}_{{ab}}^{\mathrm{Kerr}}+\eta {\chi }^{2}{h}_{{ab}},\end{eqnarray} \tag{ 71 }$$

where ${g}_{{ab}}^{\mathrm{Kerr}}$ is the Kerr metric in BL coordinates (see Equation (56)), and

$$\begin{eqnarray}{h}^{00} & = & -{(1-2M/r)}^{-1}\,[(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{1}(r)]\\ {h}^{11} & = & -(1-2M/r)[(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{1}(r)],\end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray}{h}^{22} & = & {r}^{-2}\,[(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{2}(r)]\\ {h}^{33} & = & {r}^{-2}{\sin }^{-2}\theta \,[(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{2}(r)].\end{eqnarray} \tag{ 73 }$$

Here, M is the mass of the star, $\chi =J/M$ is the dimensionless angular momentum and η parametrizes the mass quadrupole moment $q=-{\chi }^{2}(1+\eta )$, so that η = 0 corresponds to the quadrupole moment of the Kerr spacetime. The functions ${{ \mathcal F }}_{\mathrm{1,2}}$ are given in Glampedakis & Babak (2006).