Black Hole Images as Tests of General Relativity: Effects of Spacetime Geometry

We begin with the Kerr metric (Kerr 1963), which is considered the most astrophysically relevant black hole solution and describes the exterior of a static and axisymmetric black hole in General Relativity (GR). It also serves as the reference solution, to which we compare our non-GR black hole solutions. In Boyer–Lindquist (oblate spheroidal) coordinates (Boyer & Lindquist 1967), the Kerr metric line element is given by:

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -\left(1-\displaystyle \frac{2\,{r}_{{\rm{g}}}\,r}{{\rm{\Sigma }}}\right)\,{c}^{2}{{dt}}^{2}-\displaystyle \frac{4a\,{r}_{{\rm{g}}}\,r{\sin }^{2}\theta }{{\rm{\Sigma }}}\,c\,{dt}\,d\phi \\ & & +\,\displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}}\,{{dr}}^{2}+{\rm{\Sigma }}\,d{\theta }^{2}+\displaystyle \frac{{ \mathcal A }{\sin }^{2}\theta }{{\rm{\Sigma }}}\,d{\phi }^{2},\end{array}\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\rm{\Sigma }}:= {r}^{2}+{a}^{2}{\cos }^{2}\theta ,\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&{\rm{\Delta }}:= {r}^{2}-2\,{r}_{{\rm{g}}}\,r+{a}^{2},\end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray}&&{ \mathcal A }:= {\left({r}^{2}+{a}^{2}\right)}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta .\end{eqnarray} \tag{ 2c }$$

Herein, r_g ≡ GMc⁻² denotes the gravitational radius of the black hole, where M is its mass, and G and c denote Newton's gravitational constant and the speed of light, respectively. Furthermore, a ≡ J/(cM) denotes the black hole's dimensional spin parameter (units of length) and J denotes its total angular momentum. The dimensionless spin parameter may then be defined as a_* ≡ cJ/(GM²) ≡ a/ r _g. The numerical calculations in this study adopt the geometrical unit convention (wherein G = c = 1) and also let M = 1, which is equivalent to normalizing all length scales to units of r_g.

The second black hole solution we consider is the Einstein–Maxwell–Dilaton–Axion (EMDA) metric. It is chosen to serve as a particular, demonstrative non-GR black hole solution, with specific metric parameters corresponding to physical field couplings. Following García et al. (1995), the EMDA line element for a static, axisymmetric black hole may be written as:

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -\left(\displaystyle \frac{\widehat{{\rm{\Delta }}}-{a}^{2}{\sin }^{2}\theta }{\widehat{{\rm{\Sigma }}}}\right){c}^{2}{{dt}}^{2}\\ & & -\,\displaystyle \frac{2a\left(\delta -\widehat{{\rm{\Delta }}}W\right){\sin }^{2}\theta }{\widehat{{\rm{\Sigma }}}}\,c\,{dt}\,d\phi +\displaystyle \frac{\widehat{{\rm{\Sigma }}}}{\widehat{{\rm{\Delta }}}}\,{{dr}}^{2}\\ & & +\widehat{{\rm{\Sigma }}}\,d{\theta }^{2}+\displaystyle \frac{\widehat{{ \mathcal A }}{\sin }^{2}\theta }{\widehat{{\rm{\Sigma }}}}\,d{\phi }^{2},\end{array}\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&W:= 1+\left[{\beta }_{{ab}}\left(2\cos \theta -{\beta }_{{ab}}\right)+{\beta }_{a}^{2}\right]{\csc }^{2}\theta ,\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&\widehat{{\rm{\Sigma }}}:= {\rm{\Sigma }}-\left({\beta }^{2}+2{br}\right)+{r}_{{\rm{g}}}^{2}\,{\beta }_{b}\left({\beta }_{b}-2{a}_{* }\cos \theta \right),\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&\widehat{{\rm{\Delta }}}:= {\rm{\Delta }}-\left({\beta }^{2}+2{br}\right)-{r}_{{\rm{g}}}\left({r}_{{\rm{g}}}+2b\right){\beta }_{b}^{2},\end{eqnarray} \tag{ 4c }$$

$$\begin{eqnarray}&&\widehat{{ \mathcal A }}:= {\delta }^{2}-{a}^{2}\widehat{{\rm{\Delta }}}{W}^{2}{\sin }^{2}\theta ,\end{eqnarray} \tag{ 4d }$$

$$\begin{eqnarray}&&\delta := {r}^{2}-2b\,r+{a}^{2}.\end{eqnarray} \tag{ 4e }$$

Here b and β denote the coupling parameters of the dilaton and axion fields, respectively, and have units of length. For clarity, we have also defined:

$$\begin{eqnarray}&&{\beta }_{a}\equiv \displaystyle \frac{{\beta }_{* }}{{a}_{* }},\qquad {\beta }_{b}\equiv \displaystyle \frac{{\beta }_{* }}{{b}_{* }},\qquad {\beta }_{{ab}}\equiv \displaystyle \frac{{\beta }_{* }}{{a}_{* }{b}_{* }},\end{eqnarray} \tag{ 5 }$$

where b_* ≡ b/r_g and β_* ≡ β/r_g are dimensionless counterparts of these coupling parameters.

Inspecting Equation (5), one immediately notices that for the effects of the axion field to be nonzero, the dilaton coupling and the black hole spin parameter must both be nonzero, i.e., one can obtain neither a spherically symmetric nor an axisymmetric solely axion black hole solution. By contrast, through setting the axion coupling to zero, an axisymmetric dilaton black hole solution is obtained. Furthermore, since the term W is always multiplied by a, the spherically symmetric EMDA black hole solution is recovered.

The final metric we consider is the Johannsen–Psaltis (JP) metric, which is a general parameterized metric that describes the exterior solution of a rapidly spinning black hole. The JP metric introduces parametric deviations to the Kerr metric through adjustable "deviation parameters." At lowest order, four such parameters exist, of which three are of physical relevance to this study. The topological structure and geodesic integrability of the JP metric is well studied in the literature, providing a reliable platform upon which to perform observational tests of astrophysical black holes (see Johannsen 2013a; Medeiros et al. 2020). This enables the investigation of electromagnetic radiation produced in the vicinity of the event horizons of rapidly spinning black holes that cannot be described by the Kerr solution, nor be admitted as a solution of the Einstein field equations of GR. In Boyer–Lindquist–like coordinates, the JP line element (Johannsen 2013b) is written as:

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -\displaystyle \frac{\widetilde{{\rm{\Sigma }}}\,{ \mathcal B }}{{ \mathcal F }}\,{c}^{2}{{dt}}^{2}-\displaystyle \frac{2a\,\widetilde{{\rm{\Sigma }}}\,{ \mathcal C }{\sin }^{2}\theta }{{ \mathcal F }}\,c\,{dt}\,d\phi \\ & & +\,\displaystyle \frac{\widetilde{{\rm{\Sigma }}}}{{\rm{\Delta }}{A}_{5}}\,{{dr}}^{2}+\widetilde{{\rm{\Sigma }}}\,d{\theta }^{2}+\displaystyle \frac{\widetilde{{\rm{\Sigma }}}\,{ \mathcal D }{\sin }^{2}\theta }{{ \mathcal F }}\,d{\phi }^{2},\end{array}\end{eqnarray} \tag{ 6 }$$

where:

$$\begin{eqnarray}&&{ \mathcal B }:= {\rm{\Delta }}-{a}^{2}{A}_{2}^{2}{\sin }^{2}\theta ,\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{ \mathcal C }:= \left({r}^{2}+{a}^{2}\right){A}_{1}{A}_{2}-{\rm{\Delta }},\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&{ \mathcal D }:= {\left({r}^{2}+{a}^{2}\right)}^{2}{A}_{1}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta ,\end{eqnarray} \tag{ 7c }$$

$$\begin{eqnarray}&&{ \mathcal F }:= {\left[\left({r}^{2}+{a}^{2}\right){A}_{1}-{a}^{2}{A}_{2}{\sin }^{2}\theta \right]}^{2},\end{eqnarray} \tag{ 7d }$$

and the useful identity ${ \mathcal F }{\rm{\Delta }}\equiv { \mathcal B }{ \mathcal D }+{a}^{2}{{ \mathcal C }}^{2}{\sin }^{2}\theta$ holds. Finally, the terms that depend on the Kerr metric deformation parameters are defined as follows:

$$\begin{eqnarray}&&\widetilde{{\rm{\Sigma }}}:= {\rm{\Sigma }}+{r}_{{\rm{g}}}^{2}\,\sum _{n=3}^{\infty }{\epsilon }_{{\rm{n}}}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{n-2},\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&{A}_{1}:= 1+\sum _{n=3}^{\infty }{\alpha }_{1{\rm{n}}}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{n},\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&{A}_{2}:= 1+\sum _{n=2}^{\infty }{\alpha }_{2{\rm{n}}}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{n},\end{eqnarray} \tag{ 8c }$$

$$\begin{eqnarray}&&{A}_{5}:= 1+\sum _{n=2}^{\infty }{\alpha }_{5{\rm{n}}}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{n}.\end{eqnarray} \tag{ 8d }$$

The metric deformation parameters are all dimensionless and ₂ = α₁₂ = 0 in the above expressions. As noted in Johannsen (2013b), this form of the metric satisfies asymptotic flatness, recovers the correct Newtonian limit, and satisfies current PPN constraints. At lowest order in the expansion, the metric depends only on α₁₃, α₂₂, α₅₂, and ₃. When these parameters are set to zero, the Kerr metric is recovered. In this study, we let α₅₂ = 0, as it modifies only the g_rr component of the metric and fixing it to zero ensures that the (outer) event horizon radius is always equal to that of the Kerr metric.

In choosing a parameterized metric, there are several desirable conditions for black hole spacetimes to be pathology-free, namely that they are: (i) stationary, (ii) axisymmetric, (iii) satisfy asymptotic flatness, and (iv) contain no singularities or closed time-like curves external to the outer event horizon. To date, black hole solutions in alternative theories of gravity that satisfy these conditions, i.e., do not exhibit pathological behavior, generally admit a rank-2 Killing tensor and possess a fourth "Carter-like" integral of motion (Carter 1968; see Vigeland et al. 2011, for further information). The existence of a rank-2 Killing tensor is a desirable property of black hole spacetimes, ensuring the geodesic motion is not chaotic. While this restriction is not strictly necessary, deviations from geodesic separability, if present, are likely to be very small (e.g., Yagi et al. 2012) and have a negligible effect on black hole images.

The existence of a rank-2 Killing tensor ensures the geodesic motion is not chaotic (see, e.g., Zelenka & Lukes-Gerakopoulos 2017, for an example of chaotic motion in the JP spacetime). In this study we consider the subset of parameterized metrics that satisfy the integrability condition, thereby ensuring the geodesic motion remains nonchaotic. We note that this integrability condition is not strictly necessary, and mapping to all known modified black hole solutions requires this condition be relaxed (e.g., Yagi et al. 2012). The effects of nonintegrability on photon orbits have been investigated in recent studies (see, e.g., Pappas & Glampedakis 2018; Kostaros & Pappas 2022).

Several different parameterized spacetime metrics exist in the literature (Manko & Novikov 1992; Glampedakis & Babak 2006; Vigeland & Hughes 2010; Johannsen & Psaltis 2011; Vigeland et al. 2011; Rezzolla & Zhidenko 2014; Konoplya et al. 2016), all differing in the extent to which the above four conditions are satisfied. In subsequent works, some authors have extended some of these parameterized metrics to also ensure that the spacetimes yield either a separable Hamilton–Jacobi equation (e.g., the JP metric employed in this study; Johannsen 2013c) or separability in both the Hamilton–Jacobi and Klein–Gordon equations (e.g., the Konoplya–Rezzolla–Zhidenko metric; Konoplya et al. 2018; Konoplya & Zhidenko 2021). Owing to the above considerations and the motivations at the start of Section 2.3, we choose the JP metric to model parametric deviations from astrophysical Kerr black holes in this study.

We note that the covariant plasma model we have developed in Paper I and employed throughout the present work does not involve nor rely upon the existence of a fourth integral of motion. The only aspect of our calculation that is affected by such an integral of motion is the ray tracing. Several earlier studies have performed calculations of ray tracing in spacetimes that do not have a Carter-like constant. They yielded conclusions regarding the size of the black hole shadow, its relation to the size of the UPO, and its nearly circular shape that are identical to ours (see, e.g., Johannsen & Psaltis 2010a, who used the quasi-Kerr spacetime by Glampedakis & Babak 2006). The only differences introduced to the shadow shape by the nonintegrability of the spacetime are distinctive yet small deviations from the nearly circular shape in spacetimes with large deviations from Kerr (Kostaros & Pappas 2022), none of which alter the conclusions of this study.

Black hole spacetimes possess several different characteristic radii. These radii define regions that characterize different physical properties of, and processes occurring in, a given spacetime geometry. This study is concerned with three specific characteristic radii: the event horizon, r_H, the unstable photon orbit (UPO), r_UPO, and the ISCO, r_ISCO.

The covariant plasma model detailed in Section 3 specifies a radial four-velocity profile that is affected by the location of the ISCO. The radial velocity, in turn, determines the plasma number density and magnetic field strength. The critical impact parameter of geodesics comprising a given black hole image, i.e., the apparent size of the black hole shadow boundary curve, is dependent on r_UPO. When calculating geodesic motion in arbitrary spacetime geometries, a priori knowledge of the event horizon radius is necessary to specify an appropriate numerical cutoff radius, r_cut, for numerical geodesic integration algorithms. ⁵ Owing to these considerations, in this study, it is necessary to determine these radii to very high precision.

In integrating the geodesic equations of motion, we typically employ a fourth-order Runge–Kutta-Fehlberg (RKF) algorithm with adaptive step-size control. In situations where the observer is placed much farther from the black hole (fiducial distance of ∼10⁴ r_g), e.g., in the large UPO cases considered later in this paper, we employ an eighth-order RKF method with adaptive step-size control. Throughout this study, we specify an integration tolerance of 10⁻¹².

Appendix D presents the results of a cross-validation between two independent relativistic radiative transfer codes, where we show that the leading discrepancy between the codes arises from the numerical accuracy by which the ISCO radius is determined. For this reason, we calculate all three critical radii to machine double precision, i.e., ≲10⁻¹⁶, prior to geodesic integration. We detail in Appendix A the numerical procedure we use for determining the characteristic radii and present isocontours of these radii for different coupling parameters of the EMDA and JP spacetimes.

In the subsequent exploration of black hole images and their dependence on the metric properties, we will first utilize a toy emissivity model where we allow the 1.3 mm (230 GHz) emissivity to have a power-law dependence on the coordinate radius and an arbitrary scale height h/r:

$$\begin{eqnarray}&&j(r,\theta )={j}_{0}\,{r}^{-n}\exp \left\{-\displaystyle \frac{1}{2}{\left[\displaystyle \frac{\theta -\pi /2}{(h/r)\pi /2}\right]}^{2}\right\}.\end{eqnarray} \tag{ 9 }$$

This equation simplifies to a spherically symmetric emissivity as h/r goes to infinity. This will allow us to separate the effects of the spacetime from those of the additional effects introduced by the plasma model and the relativistic Doppler shifts introduced by the motion of the gas. In this model, the plasma is considered to be in freefall. Explicit expressions for the four-velocities of free-falling particle geodesics in the EMDA and JP spacetimes are presented in Appendix B.

In the bulk of the results, we will employ an analytic model that is based on the solution of the conservation laws for mass, momentum, and energy that govern the dynamics and thermodynamics of the gas accreting toward the black hole. This covariant semianalytic model is axisymmetric so that there is no dependence of any quantity on the azimuthal angle ϕ. We review here the basic equations for completeness and refer the reader to Paper I for the details and derivations.

3.2.1. Plasma Model Thermodynamic Properties

We set the equatorial electron density profile to

$$\begin{eqnarray}&&{n}_{{\rm{e}}\ ,\mathrm{eq}}(\varpi )=\displaystyle \frac{\dot{M}}{4\pi \sqrt{-g}\ (h/r)\ {m}_{{\rm{p}}}\left(-{u}^{r}\right)},\end{eqnarray} \tag{ 10 }$$

where $\varpi \equiv r\sin \theta$ is the equatorial radius, m_p is the proton mass (assuming a fully ionized hydrogen plasma), and u^r is the r − component of the plasma four-velocity. The factor related to the determinant of the metric, $g=\det {g}_{\mu \nu }$, is evaluated at the same equatorial radius. We multiply this equatorial density profile n_{e, eq} with an exponential in the polar angle θ, i.e.,

$$\begin{eqnarray}&&{n}_{{\rm{e}}}(r,\theta )={n}_{{\rm{e}},\ \mathrm{eq}}(\varpi )\exp \left\{-\displaystyle \frac{1}{2}{\left[\displaystyle \frac{\theta -\pi /2}{(h/r)\pi /2}\right]}^{m}\right\},\end{eqnarray} \tag{ 11 }$$

where the index m determines the slope of the vertical density profile. In the Newtonian limit and for an ion temperature that is constant with height, we find that m = 2 and

$$\begin{eqnarray}&&\displaystyle \frac{h}{r}=\displaystyle \frac{1}{{{ru}}^{\phi }}{\left(\displaystyle \frac{P}{\rho }\right)}^{1/2}=\sqrt{(\hat{\gamma }-1)\zeta }.\end{eqnarray} \tag{ 12 }$$

The parameters $\hat{\gamma }$ and ζ arise from energy conservation arguments and will be introduced below. We will use these expressions hereafter, unless specified otherwise.

Solving the energy conservation equation leads to the following expression for the ion temperature:

$$\begin{eqnarray}&&{T}_{{\rm{i}}}(r,\theta )=\displaystyle \frac{{m}_{p}{c}^{2}}{{k}_{{\rm{B}}}}\displaystyle \frac{R(\hat{\gamma }-1)}{(R+1)}{ \mathcal V },\end{eqnarray} \tag{ 13 }$$

where k_B is the Boltzmann constant, $\hat{\gamma }$ is the effective adiabatic index of the plasma, R is the ratio of the electron to ion temperature, and ${ \mathcal V }$ is a density-weighted integral of the dissipation function for the process that is responsible for heating the flow.

Because the flow is radiatively inefficient, the ion temperature at any radius becomes comparable to the local virial temperature, evaluated appropriately for each spacetime. Following the procedure outlined in Appendix A of Paper I and keeping only the leading-order corrections introduced by the various metric deviation parameters, we write

$$\begin{eqnarray}&&{T}_{{\rm{i}}}(r,\theta )=\displaystyle \frac{{m}_{p}{c}^{2}}{{k}_{{\rm{B}}}}\displaystyle \frac{R(\hat{\gamma }-1)}{(R+1)}\zeta \left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right){{ \mathcal T }}_{{\rm{c}}},\end{eqnarray} \tag{ 14 }$$

where ${{ \mathcal T }}_{{\rm{c}}}$ is a spacetime-dependent correction to the temperature, which is equal to unity for the Kerr black hole. For the EMDA and JP metrics used in this study, these corrections are given, respectively, by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal T }}_{{\rm{c}}} & = & 1+\displaystyle \frac{1}{3}\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)\left[7{b}_{* }+4\left({b}_{* }+1\right){\beta }_{b}^{2}\right]\\ & & +\,\displaystyle \frac{1}{9}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{2}\left[{b}_{* }\left(45{b}_{* }+51.5\right)\right.\\ & & +\,\left.\left({b}_{* }+1\right)\left(42{b}_{* }+11\right){\beta }_{b}^{2}\right]+{ \mathcal O }\left({r}^{-3}\right),\end{array}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{{ \mathcal T }}_{{\rm{c}}}=1+\displaystyle \frac{1}{6}{\left(\displaystyle \frac{{r}_{{\rm{g}}}}{r}\right)}^{2}\left(10{\alpha }_{13}-{\alpha }_{52}-5{\epsilon }_{3}\right)+{ \mathcal O }\left({r}^{-3}\right).\end{eqnarray} \tag{ 16 }$$

In both cases ζ is an order-unity factor. We then write the electron temperature as

$$\begin{eqnarray}&&{T}_{{\rm{e}}}(r,\theta )=\displaystyle \frac{{T}_{i}(r,\theta )}{R}.\end{eqnarray} \tag{ 17 }$$

Finally, we specify the magnetic field everywhere such that the plasma-β parameter is constant throughout the flow, i.e., such that

$$\begin{eqnarray}&&B(r,\theta )\propto {\left[{n}_{{\rm{e}}}(r,\theta )\ {T}_{i}(r,\theta )\right]}^{1/2}.\end{eqnarray} \tag{ 18 }$$

In the absence of synchrotron self absorption, which is negligible for the frequency and range of accretion rates of interest here, the overall normalization of the magnetic field can be specified at a fiducial equatorial location and scaled according to Equation (18). We, therefore, write

$$\begin{eqnarray}&&B\left(r,\theta \right)={B}_{0}{\left[\displaystyle \frac{{n}_{{\rm{e}}}(r,\theta )\,{T}_{{\rm{i}}}\left(r,\theta \right)}{{n}_{{\rm{e}}}\left({r}_{0},\pi /2\right)\,{T}_{{\rm{i}}}\left({r}_{0},\pi /2\right)}\right]}^{1/2},\end{eqnarray} \tag{ 19 }$$

such that B₀ is the strength of the magnetic field at the spherical radius ${r}_{0}\equiv {r}_{\mathrm{ISCO}}\left({a}_{* }=0,\theta =\pi /2\right)$, i.e., at the equatorial ISCO radius for a nonspinning black hole in the spacetime being considered. This particular normalization for non-Kerr spacetimes is chosen in order to smoothly recover the Kerr expression in Paper I, while properly accounting for cases of large metric deviation parameters where the ISCO radius for a_* = 0 can deviate significantly from the 6 r_g Schwarzschild value.

We use the analytic fitting formula for the angle-averaged emissivity derived by Mahadevan et al. (1996), which is accurate to within 2.6% for all temperatures and frequencies of interest:

$$\begin{eqnarray}&&{j}_{\nu }=\displaystyle \frac{{n}_{{\rm{e}}}\,{e}^{2}\nu }{\sqrt{3}\,c\,{K}_{2}(1/{{\rm{\Theta }}}_{e})}\,M({x}_{M}),\end{eqnarray} \tag{ 20 }$$

with M(x_M ) given by:

$$\begin{eqnarray}\begin{array}{rcl}M({x}_{M}) & = & \displaystyle \frac{4.0505\,{\mathscr{a}}}{{x}_{M}^{1/6}}\left(1+\displaystyle \frac{0.40\,{\mathscr{b}}}{{x}_{M}^{1/4}}+\displaystyle \frac{0.5316\,{\mathscr{c}}}{{x}_{M}^{1/2}}\right)\\ & & \times \,\exp \left(-1.8896\,{x}_{M}^{1/3}\right).\end{array}\end{eqnarray} \tag{ 21 }$$

Here, ${\nu }_{b}\equiv {eB}/\left(2\pi {m}_{{\rm{e}}}c\right)$ is the cyclotron frequency and

$$\begin{eqnarray}&&{x}_{M}\equiv \displaystyle \frac{2\nu }{3\,{\nu }_{b}\,{{\rm{\Theta }}}_{{\rm{e}}}^{2}},\end{eqnarray} \tag{ 22 }$$

where ${{\rm{\Theta }}}_{{\rm{e}}}\equiv {k}_{{\rm{B}}}{T}_{{\rm{e}}}/\left({m}_{{\rm{e}}}{c}^{2}\right)$ is the dimensionless electron temperature, e is the electron charge, m_e is the electron mass, and K₂(x) is the modified Bessel function of the second kind of order two and with argument x. The best-fit values of the coefficients ${\mathscr{a}},{\mathscr{b}}$, and ${\mathscr{c}}$ for different temperatures are given in Mahadevan et al. (1996).

3.2.2. Plasma Model Four-velocity Profile

For each spacetime, the equatorial radius of the ISCO provides a natural separatrix to specify the plasma velocities. On the equatorial plane and outside the ISCO, we set the azimuthal component of the four-velocity, u^ϕ , equal to the orbital velocity of test particles at the same location. The radial velocity profile depends on the efficiency of the angular momentum transport by material and magnetic stresses. In order to allow for a general form that does not depend on the specifics of angular momentum transport, we write the radial velocity as:

$$\begin{eqnarray}&&{u}_{\mathrm{eq}}^{r}(r)=-\eta {\left(\displaystyle \frac{r}{{r}_{\mathrm{ISCO}}}\right)}^{-{n}_{r}},\end{eqnarray} \tag{ 23 }$$

where η and n_r are free parameters.

On the equatorial plane and inside the ISCO, we calculate the azimuthal and radial velocity profiles by following the trajectories of test particles that freefall with the energy and angular momentum of the plasma at the ISCO radius. Throughout the flow, we assume that the polar component of the four-velocity is zero, i.e., u^θ = 0, and calculate u^t by imposing the condition g_{μ ν} u^μ u^ν = − 1.

In order to model the plasma velocities off the equatorial plane, we use the fact borne out from semianalytic models and GRMHD simulations that the azimuthal components u^ϕ are approximately constant on spherical surfaces and that the radial components u^r are approximately constant on cylindrical surfaces. We again refer the reader to Paper I for the details of the model.

We choose a set of fiducial values for the plasma parameters to employ for the majority of the paper that are representative of disk accretion and are also consistent with the numerical solutions obtained in GRMHD models of M87. In particular, unless stated otherwise, we use η = 0.1 and n_r = 1.5 for the radial velocity profile (Equation (23)), ζ = 0.25 and $\hat{\gamma }=5/3$ for the ion temperature, R = 5 for the electron temperature (Equation (17)), and B₀ = 20 G for the magnetic field scale (Equation (19)). These values lead to h/r ≃ 0.4 for the disk scale height (Equation (12)). Although we have demonstrated the lack of sensitivity of our conclusions on the specific plasma parameters in Paper I, we nevertheless use a second set of plasma parameters to show that this lack of sensitivity is not specific to the Kerr metric.

The geodesic equations of motion and the equations of general-relativistic radiative transfer (GRRT) are solved using the BHOSS code (Younsi et al. 2012, 2016, 2020). In this study, the equations of motion governing null geodesics (hereafter rays) are integrated as:

$$\begin{eqnarray}&&\displaystyle \frac{{{dk}}^{\mu }}{d\lambda }=-{{{\rm{\Gamma }}}^{\mu }}_{\alpha \beta }\,{k}^{\alpha }{k}^{\beta },\qquad \displaystyle \frac{{{dx}}^{\mu }}{d\lambda }={k}^{\mu },\end{eqnarray} \tag{ 24 }$$

where λ is the affine parameter parameterizing the ray, k^μ is the ray's four-momentum, and x^μ is its position four-vector. After specifying the expressions for the covariant metric tensor components, g_{μ ν} , the Christoffel symbols, ${{{\rm{\Gamma }}}^{\mu }}_{\alpha \beta }$, are computed via centered finite differences, with the contravariant metric tensor computed via LU-decomposition of g_{μ ν} . Equations (24) are typically integrated using a fourth-order RKF method with adaptive step-size control, and in instances where the spacetime under consideration strongly deviates from Kerr and the required precision is necessarily higher, a sixth-order RKF method is used. In numerically integrating the black hole shadow boundary curves using a bisection method (see Younsi et al. 2016), an eighth-order RKF method is used.

After determining the ray trajectories, we solve the GRRT equations along these rays, accounting for emission and the effects of attenuation of the ray intensities by the intervening media between the black hole event horizon and the observer. The equation expressing the covariant radiative transport of unpolarized, unscattered radiation (see, e.g., Lindquist 1966; Fuerst & Wu 2004; Younsi et al. 2012) is written as:

$$\begin{eqnarray}&&\displaystyle \frac{d{{ \mathcal I }}_{\nu }}{d\lambda }=-{k}^{\alpha }{u}_{\alpha }\left(-{\chi }_{\nu ,0}\,{{ \mathcal I }}_{\nu }+\displaystyle \frac{{j}_{\nu ,0}}{{{\nu }_{0}}^{3}}\right),\end{eqnarray} \tag{ 25 }$$

where u^α is the plasma four-velocity, ν denotes the observing frequency, subscript "0" denotes quantities evaluated in the comoving frame of the accretion flow, χ_ν and j_ν denote the frequency-dependent absorption and emission coefficients, respectively, and ${{ \mathcal I }}_{\nu }\equiv {I}_{\nu }/{\nu }^{3}$ is the (frequency-dependent) Lorentz-invariant intensity and I_ν the corresponding specific intensity.

We solve these equations numerically in decoupled form, as originally expressed in Younsi et al. (2012), which, for an observer (subscript "obs") at λ_obs, may be written as:

$$\begin{eqnarray}&&\displaystyle \frac{d{\tau }_{\nu }}{d\lambda }={g}^{-1}{\chi }_{\nu ,0},\end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{{ \mathcal I }}_{\nu }}{d\lambda }={g}^{-1}\left(\displaystyle \frac{{j}_{\nu ,0}}{{{\nu }_{0}}^{3}}\right){{\rm{e}}}^{-{\tau }_{\nu }},\end{eqnarray} \tag{ 26b }$$

where $g\equiv \nu /{\nu }_{0}=({k}^{\beta }{u}_{\beta }{| }_{{\lambda }_{\mathrm{obs}}})/({k}^{\alpha }{u}_{\alpha }{| }_{\lambda })$ is the relative energy shift of the photon between the observer and comoving frames, and τ_ν denotes the frequency-dependent optical depth. We solve these equations using a simple first-order Euler method, with the step size determined from the RKF geodesic integration. This formulation has the advantage of integrating the GRRT equation in tandem with the geodesic equations, i.e., from observer to source. This avoids storing geodesics in memory for the GRRT integration, enables the tracking and truncation of ray optical depth, and is computationally fast and extendable to multifrequency integration. In this study, we assume that the emission is optically thin, i.e., χ_ν = 0, and therefore integrate Equation (26b) alone.

We calculate images using the power-law emissivity profile of Equation (9) for a power-law index n = 2, a disk scale height h/r = 0.25, and an observer inclination of 45°. We also assume that the plasma is free-falling (see Appendix B for the derivation of the relevant four-velocities). Figure 1 shows horizontal image cross sections for different parameters of the EMDA and JP metrics and for different values of the black hole spin. We show the locations of the critical impact parameters as vertical dashed lines.

Zoom In Zoom Out Reset image size

Comparing these cross sections to those in Figure 1 of Paper I for the Kerr metric, we see that all of the same universal features persist for images in the non-Kerr case. In particular, the peak of the emission always occurs very close to the critical impact parameter ⁶ even when the latter varies substantially with the metric parameters. Inside the critical impact parameter, there is always a sharp brightness depression that would be identified with the black hole shadow. Finally, in the effective absence of azimuthal velocities in the free-falling regime, brightness asymmetries in the image are caused almost exclusively by frame-dragging effects. These asymmetries are relatively minor unless frame dragging is enhanced substantially beyond the Kerr value (e.g., by increasing the α₂₂ parameter in the JP metric to significantly large values).

For the remainder of this paper, we employ the full covariant plasma model to describe the accreting material and its emission characteristics around the black hole (see Section 3.2). In addition, we incorporate realistic velocities, which include substantial azimuthal components that can give rise to images with a large brightness asymmetry.

Figure 2 shows a selection of model images with the full covariant plasma model, for different values of the EMDA and JP metric parameters as well as for different black hole spins, observer inclination angles, and plasma parameters. In all cases, the images of these optically thin accretion flows are ring-like or crescent-like, with a deep central brightness depression, despite the significantly different spacetime geometries being explored. The main effect of changing the plasma parameters (see first and second rows in the figure) is to alter the width of the emission ring: shallow density profiles lead to broader rings and vice versa. Increasing the observer inclination (see the top two and bottom two rows) increases the brightness asymmetry between the approaching and receding part of the image, which is caused by relativistic Doppler effects arising from the azimuthal velocity component of the accreting material. For the model parameters displayed here, low inclinations correspond to ring-like images, whereas high inclinations correspond to crescents (see Medeiros et al. 2022 for exceptions to this behavior).

Zoom In Zoom Out Reset image size

The diameters of the bright image rings (or crescents) do not change appreciably when the black hole metric is fixed, i.e., moving down any column in this figure. On the other hand, changing the metric or its parameters, i.e., comparing different columns in the figure, alters the diameters of the bright rings. However, in all cases the diameters of the images scale with the diameters of the shadows, which are shown as dashed red lines. As a result, measuring the diameter of the ring for a black hole of known mass can indeed be used as a test of the metric (Psaltis et al. 2020).

In Figure 3, we explore in more detail the image brightness near the critical impact parameters for a number of EMDA metrics with no spin. As in the case of the analytic emissivity models discussed in Section 4.1, strong gravitational lensing near the UPOs causes a sharp increase in the brightness at impact parameters close to the critical values (denoted in the figure by vertical dashed lines). Even when the locations of the critical impact parameters change as the metric properties are varied, the dominant brightness peaks are displaced concurrently.

Zoom In Zoom Out Reset image size

Finally, we present in Figure 4 image brightness cross sections of EMDA and JP metrics with different deviation parameters, black hole spins, observer inclination angles, and plasma model parameters. This broader exploration confirms that:

• 1.  
  the diameter of the bright ring follows the shadow diameter closely,
• 2.  
  it is the critical impact parameter that plays the dominant role in determining the properties of the image, and
• 3.  
  images from optically thin accretion flows remain narrow even in non-Kerr spacetimes.

Zoom In Zoom Out Reset image size

In the following subsections, we explore a comprehensive library of non-GR black hole images, considering the variation of 10 different model parameters and subsequently quantifying the relationship between image properties and spacetime metric parameters.

Zoom In Zoom Out Reset image size

The close proximity of the peak brightness of the black hole image to the critical impact parameter makes it possible to use the characteristics of black hole images to perform measurements of the underlying spacetime properties and conduct tests of the Kerr metric. Accomplishing this in a quantitative manner requires establishing a relationship, together with its corresponding uncertainty, between the diameter ${d}_{\mathrm{im}}$ of the peak brightness, which is the observed quantity, and the diameter d_sh of the black hole shadow, which probes the metric. To this end, we introduced in previous work the calibration parameter:

$$\begin{eqnarray}&&{\alpha }_{1}\equiv \displaystyle \frac{{d}_{\mathrm{im}}}{{d}_{\mathrm{sh}}},\end{eqnarray} \tag{ 27 }$$

and used numerical and analytic models within GR to measure the range of possible values of this parameter in the Kerr spacetime (see Paper I). ⁷ Here, we extend this calculation to images in non-Kerr spacetimes.

For the majority of the spacetime parameters explored in this study, the shadow is nearly circular, and its diameter can be characterized by a single number for present purposes. Using the numerical method discussed earlier, we calculate the diameter of the shadow along different azimuthal cross sections and adopt the average value as d_sh. We then measure the image diameter ${d}_{\mathrm{im}}$ using a characterization algorithm, as before, where we filter the image at the resolution of the EHT array, calculate analytically the center of each black hole shadow, and measure over a number of azimuthal cross sections the distance of the peak emission from that center (see Paper I for additional details). We then identify the diameter of the bright emission ring as being twice the value of the median distance. Within the same algorithm, we also measure the FWHM of the shadow ring by approximating the brightness distribution along each azimuthal cross section with an asymmetric Gaussian and then calculate the median value of the widths of these Gaussians.

We carry out these measurements for a wide range of 1.3 mm images obtained using the spacetimes and plasma models described in Sections 2 and 3.

For the EMDA and JP metrics, we consider a fine-grained sampling of model parameters: (b_*, β_*) for the EMDA metric, (₃, α₁₃, α₂₂) for the JP metric, and (i, a_*, η, n_r , R, B₀, ζ) for both metrics, all of which are shown in Table 1, yielding a library of ∼2 × 10⁵ images. Each parameter affects a different aspect of the image and the spacetime. Varying a_* controls the location of the characteristic radii of the spacetime and changes physical effects such as the degree of frame-dragging, and, in the case of deviations from Kerr, the multipolar structure of the black hole. The EMDA field couplings b_* and β_* alter the spacetime geometry and are explored in Figures 10–12. Deviation parameters (₃, α₁₃, α₂₂), which are varied one at a time while the others are set to zero, control the strength of deviation of the spacetime from the Kerr metric. For instance, as the value of α₂₂ increases, the quadrupole moment of the black hole and the effects of frame-dragging are significantly enhanced beyond what is possible for even an extremal Kerr black hole. This leads to highly asymmetric shadow boundary curves as well as UPO radii very close to r_H (see Figure 5 and Appendix C for examples of this).

Table 1. Ranges of the Model Parameters Used to Construct the JP Metric and EMDA Metric Image Libraries

Model Parameter	Library Values	Parameter Description
i	15°, 30°, 45°, 60°, 75°	observer inclination angle
a*	−0.9965, −0.6509, −0.3165, 0, 0.2940, 0.5585, 0.7819, 0.9428	dimensionless spin parameter
3	−6, −4, −2, 0, 2, 4, 6, 8, 10	JP metric deviation parameter
α13	−6, −4, −2, 0, 2, 4, 6, 8, 10	JP metric deviation parameter
α22	−2, 0, 2, 4, 6, 8, 10	JP metric deviation parameter
b*	−0.4, −0.2, 0, 0.2, 0.4, 0.6, 0.8, 1	EMDA metric dilaton coupling
β*	0, 0.2, 0.4, 0.6, 0.8, 1	EMDA metric axion coupling
η	0.05, 0.1, 0.2	radial velocity profile scale
nr	0.5, 1.0, 1.5	radial velocity profile index
R	5, 10	ion-to-electron temperature ratio
B0 [G]	5, 20, 50	magnetic field scale
ζ	0.25, 0.4	proxy for disk scale height a

Notes. Every image has a field of view of [−20 r_g, 20 r_g] in both the horizontal and vertical directions and is calculated at a resolution of 128 × 128 pixels. All images are produced at an observing wavelength of 1.3 mm. The total number of models in each metric library is ∼10⁵. Models with deviation parameters outside the allowed range are omitted (see Appendix. A). JP metric models fix two deviation parameters to zero while varying the third parameter. EMDA metric models vary b_* and β_* simultaneously.

^a Where in hydrostatic equilibrium and far from r_UPO the disk scale height is $h/r=\sqrt{\left(\hat{\gamma }-1\right)\zeta }$, i.e., Equation (12).

Download table as:  ASCIITypeset image

Adjusting (η, n_r ) directly alters the plasma radial four-velocity and its profile steepness as a function of radius, affecting the amount of extended emission produced in black hole images within the ISCO radius in particular, but also in its outer vicinity. The parameters R and B₀ control the radiative properties of the plasma. As noted in Paper I, we do not consider values of R = 1 since this is inconsistent with the assumption of a radiatively inefficient accretion flow, as anticipated for EHT target sources (see, e.g., EHT Collaboration 2019b). We also avoid consideration of models with R > 10, because larger values of R strongly suppress emission from the accretion flow in the models. We choose values of B₀, which, based on isothermal one-zone modeling, yield electron temperatures consistent with the observed maximum brightness temperature of M87 (EHT Collaboration 2021). The largest value of B₀ = 50 G, which is slightly larger than the upper-bound of the one-zone model, allows us to probe the broader range of potential electron number densities possible within non-Kerr accretion flow models. Finally, in all models, we fix the adiabatic index as $\hat{\gamma }=5/3$ while varying ζ, in effect adjusting the disk scale height. We find that library image morphologies exhibit the weakest dependence on ζ.

In Figure 6 we plot the fractional width for each model image against the fractional diameter difference α₁ −1. As in the case of the Kerr metric discussed in Paper I, the fractional diameter difference is almost always positive and small, i.e., the bright ring has a slightly larger diameter than the shadow in the vast majority of cases (see Figure 10 of Paper I, but note that the distributions there appear narrower because the panels are separated into positive and negative black hole spins). This is true for all metric deformation parameters in the various metrics we considered as well as for the details of the plasma model. This is a direct consequence of strong gravitational lensing effects close to r_UPO that cause a rapid increase in image brightness near the critical impact parameter and a strong brightness depression interior to it.

Zoom In Zoom Out Reset image size

The red filled regions in Figure 6 correspond to images in which the fractional width of the ring is smaller than the fractional diameter difference, while the diagonal dashed lines correspond to half of this value. In other words, models present in that region correspond to those images where the bright ring is displaced from the black hole shadow by more than the ring width. The fact that only a very small fraction of model images lie in that region shows that, even in non-Kerr metrics, ring-like images scale with the diameter of the black hole shadow and, in the vast majority of cases, are not disjoint from it.

In the case of the non-Kerr images used in calculating the α-calibration plots in Figure 6, we have restricted the analysis to images satisfying two particular requirements.

The first requirement is that a given image must contain at least half of a ring-like feature so that a well-defined diameter is measurable. As in Paper I, where we defined the fractional coverage of a circle in an image, ${ \mathcal F }$, we consider only images with ${ \mathcal F }\geqslant 0.5$.

The second requirement is that the aspect ratio of the non-Kerr black hole shadow, i.e., the ratio of its principal axes, is $0.8\leqslant {{ \mathcal R }}_{\mathrm{non}-\mathrm{Kerr}}\leqslant 1.38$. We note that for the Kerr spacetime, this aspect ratio is always constrained to $1\lesssim {{ \mathcal R }}_{\mathrm{Kerr}}\lesssim 1.15;$ our requirement corresponds to aspect ratios within 20% of Kerr.

These two requirements are based on the observation that the 1.3 mm image of the M87 black hole is a nearly circular ring (EHT Collaboration 2019c) and are imposed to eliminate image morphologies that prevent us from using a single radius to characterize images, as is required in the α-calibration. Examples of images eliminated based on the above considerations are presented and discussed in Appendix C.

In the choice of values of metric deviation parameters considered in Sections 4.3 and 4.4, we explored those parameters limited to the ranges shown in Figures 10–13. In particular, for the JP metric, we have restricted ourselves to black holes with UPO radii up to ∼5 r_g. However, in order to demonstrate the dominant role of gravitational lensing and, in particular, of the critical impact parameter in determining the characteristic features of optically thin black hole images, Figures 7 and 8 explore the images for a number of spacetime parameters with extreme deviations from Kerr. (Note that such deviations are already excluded by the present M87 images).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In particular, these figures show 1.3 mm images of black holes described by the JP metric with parameters chosen such that the radius of the equatorial circular UPO is equal to 10 r_g, 100 r_g, and 1000 r_g. For these three values of r_UPO, we choose the dimensionless spin of each of these three black holes to be −0.9965, 0, and 0.9987, yielding nine images per observer inclination angle and probing a broad range of a_* and r_UPO. For comparison, for a Kerr black hole, 1 r_g < r_UPO < 4 r_g, with r_UPO = 3 r_g in the case of a Schwarzschild black hole (a_* = 0). For each pair (a_*, r_UPO), the corresponding value of α₁₃ is obtained numerically to an accuracy better than 10⁻¹⁶. Table 2 presents the numerical values of the parameters that produce these extreme UPO radii. We note that while these values of r_UPO are unphysical in that they are much larger than what is constrained by EHT measurements (EHT Collaboration 2019c), the underlying parameter choices for the metric do not give rise to any pathological behavior or violate any fundamental properties of the spacetime manifold (e.g., closed time-like curves, curvature singularities outside the event horizon, or positive-definite metric determinant values).

Table 2. Values of the JP Deviation Parameter α₁₃ That Yield Black Hole UPO Radii of 10 r_g, 100 r_g, and 1000 r_g

${r}_{\mathrm{UPO}}\left({r}_{{\rm{g}}}\right)$	a*	α13	${r}_{\mathrm{ISCO}}\left({r}_{{\rm{g}}}\right)$
	−0.9965	3.27447 × 102	3.59824 × 101
10	0	4.11765 × 102	3.84384 × 101
	0.9987	4.90103 × 102	4.06180 × 101
	−0.9965	4.84779 × 105	1.20952 × 103
100	0	4.92386 × 105	1.21871 × 103
	0.9987	4.99890 × 105	1.22773 × 103
	−0.9965	4.98500 × 108	3.86751 × 104
1000	0	4.99249 × 108	3.87041 × 104
	0.9987	4.99998 × 108	3.87331 × 104

Note. For every UPO radius, a_* is chosen as (−0.9965, 0, 0.9987). For each calculated value of α₁₃, the resulting ISCO radius is also shown. By construction, the JP spacetime event horizon radius is the corresponding Kerr value, i.e., r_H/r_g ≃ (1.08359, 2, 1.05097).

Download table as:  ASCIITypeset image

The plasma parameters for these extreme examples are: (η, n_r , R, B₀, ζ) = (0.5, 1.5, 5, 20 G, 0.25). Each row in Figures 7 and 8 presents three images for a fixed value of a_*, varying the value of r_UPO, along with their corresponding normalized horizontal and vertical intensity cross sections. In order to demonstrate the scaling of the image diameters with the spacetime properties, both the images and the cross sections are displayed with the impact parameters divided by the corresponding radii of the UPOs. Despite such extreme deviations from the Kerr metric, which yield enormous UPO radii, one sees that the peak flux in the image is still effectively coincident with the critical impact parameter, and a characteristic pronounced central brightness depression is still seen, demonstrating that such image features are effectively independent of the accretion plasma and sensitive to only one property of the spacetime geometry, namely the location of r_UPO.

As described in Johannsen (2013b), a sufficiently general definition of the event horizon of a general black hole may be encompassed by the condition r = H(θ), which is defined by the condition:

$$\begin{eqnarray}&&{g}^{{rr}}(r,\theta )-2{g}^{r\theta }(r,\theta )\left[\displaystyle \frac{{dH}(\theta )}{d\theta }\right]+{g}^{\theta \theta }(r,\theta ){\left[\displaystyle \frac{{dH}(\theta )}{d\theta }\right]}^{2}=0,\end{eqnarray} \tag{ A1 }$$

where g^{r θ} = 0 for all metrics investigated in this study. Equation (A1) is solved numerically through discretization of the zenith derivative using a backward finite-difference approximation, starting on the polar axis (θ = 0) and iterating until θ = π/2.

For the Kerr black hole, the (outer) event horizon radius is given by the larger solution to Δ = 0, namely:

$$\begin{eqnarray}&&{r}_{{\rm{H}},\mathrm{Kerr}}={r}_{{\rm{g}}}\left(1+\sqrt{1-{a}_{* }^{2}}\right).\end{eqnarray} \tag{ A2 }$$

The EMDA black hole event horizon radius is given by:

$$\begin{eqnarray}&&{r}_{{\rm{H}},\mathrm{EMDA}}={r}_{{\rm{g}}}\left[\left(1+{b}_{* }\right)+\sqrt{\left(1+{\beta }_{b}^{2}\right){\left(1+{b}_{* }\right)}^{2}-{a}_{* }^{2}}\right].\end{eqnarray} \tag{ A3 }$$

Given a choice of metric parameters that yield a "well-behaved" metric, when α₅₂ = 0, the event horizon of the JP metric is precisely the Kerr event horizon. In this study, the event horizon radii for the EMDA and JP metrics are both calculated fully numerically via Equation (A1), with the analytic expressions above serving as a check of the accuracy of the numerical root finding algorithms used to determine r_H.

Consider a general static and axisymmetric metric:

$$\begin{eqnarray}&&{{ds}}^{2}={g}_{{tt}}{{dt}}^{2}+2{g}_{t\phi }{dt}\,d\phi +{g}_{{rr}}{{dr}}^{2}+{g}_{\theta \theta }d{\theta }^{2}+{g}_{\phi \phi }d{\phi }^{2}.\end{eqnarray} \tag{ A4 }$$

The Lagrangian equations of motion $2{ \mathcal L }={g}_{\mu \nu }{\dot{x}}^{\mu }{\dot{x}}^{\nu }$ for this metric yield the constants ${p}_{t}\equiv -E={g}_{{tt}}\dot{t}+{g}_{t\phi }\dot{\phi }$ and ${p}_{\phi }\equiv {L}_{{\rm{z}}}={g}_{t\phi }\dot{t}+{g}_{\phi \phi }\dot{\phi }$, where an overdot denotes differentiation with respect to proper time (τ). For a fluid particle, the normalization condition u_μ u^μ = −1, together with the assumption of circular orbits in the equatorial plane, i.e., θ = π/2 and $\dot{\theta }=0$, when paired with the radial equation of geodesic motion yields:

$$\begin{eqnarray}&&{\rho }^{2}\dot{t}=\left({g}_{t\phi }{L}_{{\rm{z}}}+{g}_{\phi \phi }E\right),\end{eqnarray} \tag{ A5a }$$

$$\begin{eqnarray}&&{\rho }^{2}\dot{\phi }=-\left({g}_{{tt}}{L}_{{\rm{z}}}+{g}_{t\phi }E\right),\end{eqnarray} \tag{ A5b }$$

$$\begin{eqnarray}&&{g}_{{rr}}{\dot{r}}^{2}=-\left({g}_{{tt}}+2{g}_{t\phi }{\rm{\Omega }}+{g}_{\phi \phi }{{\rm{\Omega }}}^{2}\right){\dot{t}}^{2}-1,\end{eqnarray} \tag{ A5c }$$

where ${\rm{\Omega }}:= \dot{\phi }/\dot{t}$ and ${\rho }^{2}\equiv {({g}_{t\phi })}^{2}-{g}_{{tt}}\,{g}_{\phi \phi }$.

Circular equatorial orbits imply $\dot{r}=\ddot{r}=0$, from which one obtains the condition g_tt,r + 2g_{t ϕ,r} Ω + g_{ϕ ϕ,r} Ω² = 0, where f_,μ ≡ ∂_μ f:=∂f/∂x^μ . Solving this condition yields the orbital angular velocity of circular orbits in terms of first derivatives of the metric:

$$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{-{g}_{t\phi ,r}\pm \sqrt{{({g}_{t\phi ,r})}^{2}-{g}_{{tt},r}\,{g}_{\phi \phi ,r}}}{{g}_{\phi \phi ,r}},\end{eqnarray} \tag{ A6 }$$

where the positive sign denotes orbits which are co-rotating with a prograde spinning black hole (a_* > 0). Counter-rotating orbits (negative sign in Equation (A6)) are not considered separately in this study, as we allow for the black hole spin to be negative. The condition $\dot{r}=0$ yields:

$$\begin{eqnarray}&&{u}^{t}={\left(-{g}_{{tt}}-2{g}_{t\phi }{\rm{\Omega }}-{g}_{\phi \phi }{{\rm{\Omega }}}^{2}\right)}^{-1/2},\end{eqnarray} \tag{ A7 }$$

from which the energy and angular momentum of a particle in circular orbit in the equatorial plane may be written as:

$$\begin{eqnarray}&&E=-{u}^{t}\left({g}_{{tt}}+{g}_{t\phi }{\rm{\Omega }}\right),\end{eqnarray} \tag{ A8a }$$

$$\begin{eqnarray}&&{L}_{z}={u}^{t}\left({g}_{t\phi }+{g}_{\phi \phi }{\rm{\Omega }}\right).\end{eqnarray} \tag{ A8b }$$

The UPO radius is given by the smallest real value of r ≥ r_H for which, in the limit r → r_UPO, both E⁻¹ and ${L}_{{\rm{z}}}^{-1}$ tend to zero. Consequently, r_UPO is obtained from ${({u}^{t})}^{-1}=0$ and is calculated numerically as the solution of the equation:

$$\begin{eqnarray}&&\left({g}_{{tt}}+2{g}_{t\phi }{\rm{\Omega }}+{g}_{\phi \phi }{{\rm{\Omega }}}^{2}\right){| }_{r={r}_{\mathrm{UPO}}}=0.\end{eqnarray} \tag{ A9 }$$

Substituting Equations (A5a) and (A5b) into Equation (A5c) yields, upon simplification:

$$\begin{eqnarray}&&{\rho }^{2}{g}_{{rr}}{\dot{r}}^{2}={g}_{\phi \phi }{E}^{2}+2{g}_{t\phi }{{EL}}_{{\rm{z}}}+{g}_{{tt}}{L}_{{\rm{z}}}^{2}-{\rho }^{2}.\end{eqnarray} \tag{ A10 }$$

Recalling the condition $\dot{r}=\ddot{r}=0$, upon application to (A10), one obtains:

$$\begin{eqnarray}&&{g}_{\phi \phi ,{rr}}{E}^{2}+2{g}_{t\phi ,{rr}}{{EL}}_{{\rm{z}}}+{g}_{{tt},{rr}}{L}_{{\rm{z}}}^{2}-{({\rho }^{2})}_{,{rr}}=0.\end{eqnarray} \tag{ A11 }$$

The ISCO radius is then obtained numerically as the smallest real solution of (A11) satisfying r_ISCO ≥ r_UPO.

This subsection presents figures of isocontours of UPO and ISCO radii for the JP and EMDA spacetimes, as a function of different metric parameters and a_*. The excluded regions for the JP spacetime are given by: ${\epsilon }_{3}\leqslant -{r}_{+}^{3}$, ${\alpha }_{13}\leqslant -{r}_{+}^{3}$, and ${\alpha }_{22}\leqslant -{r}_{+}^{2}$, where ${r}_{+}\equiv 1+\sqrt{1-{a}_{* }^{2}}$. For the EMDA spacetime, an exclusion region exists when β_* = 0, and is given by: ${b}_{* }\leqslant -(1-{a}_{* }^{2})/2$.

In the case of the JP metric, the parameter α₅₂ affects only the g_rr metric component, which does not affect the UPO and ISCO radii; thus, fixing α₅₂ = 0 ensures the JP and Kerr horizon radii are always coincident. The ₃ parameter only affects, albeit weakly, the ISCO radius (see Figure 10). JP metric isocontours of the ISCO and UPO radii, for varying values of one deviation parameter and spin (with all other parameters set to zero), are shown in Figures 10 and 11. These results are in excellent agreement with Johannsen (2013a), where we highlight the shaded red region in the α₁₃ UPO plot of Figure 11, which is absent in Figure 2 of Johannsen (2013a) since they considered a truncated Taylor expansion approximation of the underlying equations rather than the full numerical solution employed here. This was subsequently rectified in Figure 6 of Johannsen (2013b). For the EMDA metric, the cases β_* = 0, b_* = 1, and b_* = 0.1 are shown in Figures 12 and 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is a straightforward albeit lengthy exercise to determine the EMDA metric equations of motion from the Hamilton–Jacobi equation. After specifying an appropriate ansatz for the separation constant of the Hamilton–Jacobi equation, simplification of the resulting equations of motion of a test particle in the EMDA spacetime yields:

$$\begin{eqnarray}&&\mu \widehat{{\rm{\Sigma }}}\left(\displaystyle \frac{{dt}}{d\tau }\right)=-{aW}\left({aEW}{\sin }^{2}\theta -{L}_{{\rm{z}}}\right)+\displaystyle \frac{\delta \,{ \mathcal P }}{\widehat{{\rm{\Delta }}}},\end{eqnarray} \tag{ B1a }$$

$$\begin{eqnarray}&&\mu \displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}\left(\displaystyle \frac{{dr}}{d\tau }\right)=\pm \sqrt{\widehat{{ \mathcal R }}(r)},\end{eqnarray} \tag{ B1b }$$

$$\begin{eqnarray}&&\mu \displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}\left(\displaystyle \frac{d\theta }{d\tau }\right)=\pm \sqrt{\widehat{{\rm{\Theta }}}\left(\theta \right)},\end{eqnarray} \tag{ B1c }$$

$$\begin{eqnarray}&&\mu \displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}\left(\displaystyle \frac{d\phi }{d\tau }\right)=-\left({aEW}-{L}_{{\rm{z}}}^{2}{\csc }^{2}\theta \right)+\displaystyle \frac{a\,{ \mathcal P }}{\widehat{{\rm{\Delta }}}},\end{eqnarray} \tag{ B1d }$$

and where

$$\begin{eqnarray}\begin{array}{rcl}\widehat{{ \mathcal R }}(r) & := & {{ \mathcal P }}^{2}-\widehat{{\rm{\Delta }}}\left[{\mu }^{2}\left({r}^{2}-2{br}-{\beta }^{2}\right)+{ \mathcal Q }\right],\\ \widehat{{\rm{\Theta }}}\left(\theta \right) & := & -{\mu }^{2}\left[{a}^{2}{\cos }^{2}\theta +{r}_{{\rm{g}}}^{2}{\beta }_{b}\left({\beta }_{b}-2{a}_{* }\cos \theta \right)\right]\end{array}\end{eqnarray} \tag{ B2a }$$

$$\begin{eqnarray}&&-\displaystyle \frac{{\left({L}_{{\rm{z}}}-{aEW}{\sin }^{2}\theta \right)}^{2}}{{\sin }^{2}\theta }+{ \mathcal Q },\end{eqnarray} \tag{ B2b }$$

$$\begin{eqnarray}&&{ \mathcal P }:= E\delta -{{aL}}_{{\rm{z}}},\end{eqnarray} \tag{ B2c }$$

$$\begin{eqnarray}&&{ \mathcal Q }:= Q+{\left({L}_{{\rm{z}}}-{aE}\right)}^{2}.\end{eqnarray} \tag{ B2d }$$

The velocities of radially free-falling test particles are obtained when the constants of motion are E = μ, L_z = 0, and Q = 0. Upon substitution into Equation (B1d), one obtains:

$$\begin{eqnarray}&&\widehat{{\rm{\Sigma }}}{\left(\displaystyle \frac{{dt}}{d\tau }\right)}_{\mathrm{ff}}=-{a}^{2}{W}^{2}{\sin }^{2}\theta +\displaystyle \frac{{\delta }^{2}}{\widehat{{\rm{\Delta }}}},\end{eqnarray} \tag{ B3a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{{dr}}{d\tau }\right)}_{\mathrm{ff}}=-\sqrt{{\delta }^{2}-\widehat{{\rm{\Delta }}}\left(\delta -{\beta }^{2}\right)},\end{eqnarray} \tag{ B3b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{d\theta }{d\tau }\right)}_{\mathrm{ff}}=\pm \sqrt{-\beta \left[\beta +\displaystyle \frac{X}{{b}^{2}}+\displaystyle \frac{\beta {X}^{2}}{{a}^{2}{b}^{4}{\sin }^{2}\theta }\right]},\end{eqnarray} \tag{ B3c }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widehat{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{d\phi }{d\tau }\right)}_{\mathrm{ff}}=a\left(-W+\displaystyle \frac{\delta }{\widehat{{\rm{\Delta }}}}\right),\end{eqnarray} \tag{ B3d }$$

where $X\equiv 2\,a\,b\,{r}_{{\rm{g}}}\cos \theta +\beta ({b}^{2}-{r}_{{\rm{g}}}^{2})$ and the subscript "ff" denotes freefall. In the above, the negative root of the radial motion is taken, corresponding to freefall onto the black hole. Note, however, that even in the freefall case, the zenith motion does not vanish for the full EMDA metric, due to the axion field coupling parameter, i.e., the aforementioned ansatz fails since the axion field coupling breaks the integrability of the motion. In this study we therefore set β_* = 0 when considering freefall motion in the EMDA spacetime. This ensures Equations (B1c) and (B3c) vanish when L_z = Q = 0, thereby constraining the motion of a test particle to be confined to a two-dimensional plane, as required. The standard Kerr expressions are recovered in the limit b → 0.

Similar to the previous subsection for the EMDA spacetime, the equations of motion for particles in the JP metric may be re-expressed in the free-falling regime as:

$$\begin{eqnarray}&&\widetilde{{\rm{\Sigma }}}{\left(\displaystyle \frac{{dt}}{d\tau }\right)}_{\mathrm{ff}}=-{a}^{2}{\sin }^{2}\theta +\displaystyle \frac{{\left({r}^{2}+{a}^{2}\right)}^{2}{A}_{1}^{2}}{{\rm{\Delta }}},\end{eqnarray} \tag{ B4a }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widetilde{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{{dr}}{d\tau }\right)}_{\mathrm{ff}}=-\sqrt{{A}_{5}}\ \sqrt{{ \mathcal K }+2\,{r}_{{\rm{g}}}r\left({r}^{2}+{a}^{2}\right)},\end{eqnarray} \tag{ B4b }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widetilde{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{d\theta }{d\tau }\right)}_{\mathrm{ff}}=0,\end{eqnarray} \tag{ B4c }$$

$$\begin{eqnarray}&&\displaystyle \frac{\widetilde{{\rm{\Sigma }}}}{c}{\left(\displaystyle \frac{d\phi }{d\tau }\right)}_{\mathrm{ff}}=a\left[-1+\displaystyle \frac{\left({r}^{2}+{a}^{2}\right){A}_{1}{A}_{2}}{{\rm{\Delta }}}\right],\end{eqnarray} \tag{ B4d }$$

where ${ \mathcal K }\equiv \left({A}_{1}^{2}-1\right){\left({r}^{2}+{a}^{2}\right)}^{2}-{\rm{\Delta }}\left(\widetilde{{\rm{\Sigma }}}-{\rm{\Sigma }}\right)$. Note that for the JP metric the existence of three constants of motion is guaranteed by construction. The standard Kerr metric expressions are recovered when all deviation parameters are zero.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size