A Parametric Model for the Shapes of Black Hole Shadows in Non-Kerr Spacetimes

Much of our current understanding of black holes relies on the assumption that they are described by the Kerr solution to the Einstein equations. One possible avenue for conducting a test of this assumption with observations in the electromagnetic spectrum is by measuring the size and shape of the shadow a black hole casts on the surrounding emission (Johannsen & Psaltis 2010b; Psaltis 2018).

The black hole shadow depends only on the geometry of the spacetime and not on the astrophysics of the accretion process. For a Kerr black hole, the shadow has a radius of 5 ± 0.2 GMc⁻² for all spins and observer inclinations (see, e.g., de Vries 2000; Johannsen & Psaltis 2010b; Chan et al. 2013, see also Figure 1). Therefore, measuring the size of the shadow of a black hole of known mass constitutes a null hypothesis test of the Kerr metric (Psaltis et al. 2015). At the same time, the shape of the shadow for a Kerr black hole is nearly circular for all but the highest spins. For a general spacetime, the shape of the shadow depends primarily on the deviation of the spacetime quadrupole from its Kerr value (Johannsen & Psaltis 2010b). As a result, measuring the shadow shape leads to a test of the general relativistic no-hair theorem (Psaltis et al. 2016).

The Event Horizon Telescope (EHT) is a millimeter very long baseline interferometry experiment that has produced the first ever image of the black hole in M87 resolved at horizon scales and has measured the size and shape of its shadow (Event Horizon Telescope Collaboration et al. 2019a). The initial analysis of the EHT data has revealed no striking deviations from the predictions for the Kerr metric (Event Horizon Telescope Collaboration et al. 2019b); improved measurements with more interferometric stations and new algorithms will tighten these constraints in the near future. The second primary target of the EHT is the Galactic Center black hole, Sgr A*. This has the largest angular size on the sky of any currently known black hole (e.g., Johannsen et al. 2012) and well constrained mass and distance (Ghez et al. 2008; Gillessen et al. 2009; Gravity Collaboration et al. 2018); therefore, measuring the shape and size of its shadow is expected to provide a precise test of the Kerr metric.

In order to use the shape of the black hole shadow to test the Kerr metric with both Sgr A* and the black hole in M87, we need to explore the shapes and sizes of shadows in other metrics. Significant progress has been made in simulating the observational appearance of black holes and more exotic compact objects with different spacetimes and in different theories of gravity (see, e.g., Bambi & Freese 2009; Amarilla et al. 2010; Bambi & Yoshida 2010; Papnoi et al. 2014; Mizuno et al. 2018; Shaikh et al. 2019). Additionally, several parameterized metrics have been developed that can be used to explore the range of non-Kerr spacetimes in a way that is agnostic to the underlying physical theory (for a few axisymmetric examples, see Manko & Novikov 1992; Glampedakis & Babak 2006; Vigeland & Hughes 2010; Johannsen & Psaltis 2011; Vigeland et al. 2011; Konoplya et al. 2016). The black hole shadows that result from some of these parameterized metrics have also been explored (see, e.g., Johannsen & Psaltis 2010b; Johannsen 2013a; Ghasemi-Nodehi et al. 2015; Younsi et al. 2016; Wang et al. 2017, and Cunha & Herdeiro 2018 for a review).

Deriving a shadow from a given metric is relatively straightforward. However, using a shadow to place constraints on the metrics that could have created it presents a significant challenge as the mapping between the two is highly nonlinear. Abdujabbarov et al. (2015) proposed an algorithm that represents the shadow shapes empirically as a sum of Legendre polynomials, with no reference to the underlying metrics. Here, we use Principal Component Analysis to show that the shapes of all black hole shadows generated by several of the parameterized metrics mentioned above can be represented by a small set of functions. This allows us to generate a general parametric model for shadow shapes that is largely agnostic of the underlying spacetime metric, is more compact than a general polynomial expansion, and can be used to derive metric constraints from the EHT data.

The paper is organized as follows. In Section 2, we introduce the parameterized metrics we consider; we outline our ray-tracing simulations in Section 3. In Section 4, we introduce principal components analysis and apply it to the simulated set of shadows. In Section 5, we assess the accuracy of the PCA reconstructions and explore the relationship between the amplitude of PCA components and the metric parameters. Finally, we discuss the implications of our work in Section 6.

In this section, we introduce a number of parameterized metrics and write them explicitly in the form we use in the simulations. We select a few metrics from those that have been published, while prioritizing metrics that differ significantly from each other. For example, since the metric of Vigeland et al. (2011) builds upon the metric of Collins & Hughes (2004) and the metric of Cardoso et al. (2014) builds upon the metric of Johannsen & Psaltis (2011), we choose to only include one metric from each of these pairs. We note that many of the proposed parameterized metrics have pathologies such as naked singularities, closed timelike loops, or non-Lorentzian geometries, that can significantly complicate our numerical calculations. Resolving such pathologies is outside the scope of the present work, so we restrict ourselves to three representative metrics that have been investigated in detail for pathologies in Johannsen (2013a; in particular, see Table 1 in this reference for a summary of such pathologies).

Our starting point is the Kerr metric, which in Boyer–Lindquist coordinates takes the form (see, e.g., Bardeen et al. 1972)

$$\begin{eqnarray}\begin{array}{rcl}{g}_{{tt}}^{{\rm{K}}} & = & -\left(1-\displaystyle \frac{2{Mr}}{{\rm{\Sigma }}}\right),\\ {g}_{{rr}}^{{\rm{K}}} & = & \displaystyle \frac{{\rm{\Sigma }}}{{\rm{\Delta }}},\,\,\,{g}_{\theta \theta }^{{\rm{K}}}={\rm{\Sigma }},\\ {g}_{\phi \phi }^{{\rm{K}}} & = & \left({r}^{2}+{a}^{2}+\displaystyle \frac{2{{Ma}}^{2}r{\sin }^{2}\theta }{{\rm{\Sigma }}}\right){\sin }^{2}\theta ,\\ {g}_{t\phi }^{{\rm{K}}} & = & -\displaystyle \frac{2{Mar}{\sin }^{2}\theta }{{\rm{\Sigma }}},\end{array}\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}\equiv {a}^{2}-2{Mr}+{r}^{2},\,\,\,{\rm{\Sigma }}\equiv {r}^{2}+{a}^{2}{\cos }^{2}\theta ,\end{eqnarray} \tag{ 2 }$$

M is the black hole mass, and a is the dimensional spin parameter, a ≡ J/M. We will also make use of the dimensionless spin parameter a_* = J/M². Here and throughout the paper, we use gravitational units and set G = c = 1. In Figure 1, we show, as a point of comparison to the following figures, the effect of changing the black hole spin and the inclination of the observer on the black hole shadow of the Kerr metric. As discussed earlier, the size and shape of a black hole shadow depends very weakly on spin; Kerr shadows are approximately circular, except for extremely high values of spin.

Zoom In Zoom Out Reset image size

2.1. The Quasi-Kerr Metric

The quasi-Kerr (QK) metric is based on the work of Hartle and Thorne (Hartle 1967; Hartle & Thorne 1968) for slowly spinning neutron stars and was adapted to describe general vacuum spacetimes by Glampedakis & Babak (2006). This metric has all but the quadrupole moments equal to those of the Kerr metric; the quadrupole moment is set to

$$\begin{eqnarray}&&{Q}_{\mathrm{QK}}=-M({a}^{2}+{\epsilon }_{\mathrm{QK}}{M}^{2}),\end{eqnarray} \tag{ 3 }$$

where _QK measures the strength of the deviation from the Kerr quadrupole. The metric is perturbed in contravariant form such that

$$\begin{eqnarray}&&{g}_{\mathrm{QK}}^{\mu \nu }={g}_{{\rm{K}}}^{\mu \nu }+{\epsilon }_{\mathrm{QK}}{h}_{\mathrm{QK}}^{\mu \nu },\end{eqnarray} \tag{ 4 }$$

where ${g}_{{\rm{K}}}^{\mu \nu }$ are the metric components for the Kerr metric in contravariant form and ${h}_{\mathrm{QK}}^{\mu \nu }$ are the metric perturbations for the QK metric. In Boyer–Lindquist-like coordinates, the elements of the perturbed metric are given by (Glampedakis & Babak 2006; Johannsen 2013b)

$$\begin{eqnarray}\begin{array}{rcl}{h}_{\mathrm{QK}}^{{tt}} & = & \displaystyle \frac{r}{r-2M}(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{1}(r),\\ {h}_{\mathrm{QK}}^{{rr}} & = & \displaystyle \frac{r-2M}{r}(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{1}(r),\\ {h}_{\mathrm{QK}}^{\theta \theta } & = & -\displaystyle \frac{1}{{r}^{2}}(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{2}(r),\\ {h}_{\mathrm{QK}}^{\phi \phi } & = & -\displaystyle \frac{1}{{r}^{2}{\sin }^{2}\theta }(1-3{\cos }^{2}\theta ){{ \mathcal F }}_{2}(r),\end{array}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{1}(r) & = & -\displaystyle \frac{5(r-M)}{8{Mr}(r-2M)}(2{M}^{2}+6{Mr}-3{r}^{2})\\ & & -\,\displaystyle \frac{15r(r-2M)}{16{M}^{2}}\mathrm{ln}\left(\displaystyle \frac{r}{r-2M}\right),\\ {{ \mathcal F }}_{2}(r) & = & \displaystyle \frac{5}{8{Mr}}(2{M}^{2}-3{Mr}-3{r}^{2})\\ & & +\,\displaystyle \frac{15}{16{M}^{2}}({r}^{2}-2{M}^{2})\mathrm{ln}\left(\displaystyle \frac{r}{r-2M}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

The QK metric is a solution to the vacuum Einstein equations for small spins, a_* ≪ 1, and reduces to the Kerr metric when _QK = 0. The validity of the QK metric is limited to regions where r ≳ 2M since the logarithm in Equation (6) diverges at r = 2M. In general, this metric describes a naked singularity (Johannsen 2013a).

Johannsen & Psaltis (2010a) and Johannsen (2013b) found that, to first order in the deviation parameter, the conditions a_* ≲ 0.4 and −0.5 ≲ _QK ≲ 0.5 are necessary but not sufficient to ensure that Lorentz violations and closed timelike loops are contained within the circular photon orbit. In this work, we solve the geodesics for the full metric shown in Equations (4) and (5) without any additional approximations or expansions. For this reason, the range of allowed values for _QK is significantly smaller than the range explored in Johannsen & Psaltis (2010b), who expanded the geodesic equations to first order in the perturbation parameter. In order to avoid pathologies in our work, we require, e.g., _QK ≲ 0.2 for a black hole with i = 90° and a_* = 0.4 and _QK ≲ 0.35 for a black hole with i = 90° and a_* = 0.1.⁴

2.2. The Modified Kerr Metric

This metric proposed by Johannsen & Psaltis (2011) introduces polynomial perturbations to both the ${rr}$- and $\theta \theta$-components of the Schwarzschild metric and, following Vigeland et al. (2011), uses the Newman–Janis algorithm (Newman & Janis 1965; Drake & Szekeres 2000) to transform this into that of a rotating compact object. The resulting metric is not Ricci flat. Here we will follow the expansions of Johannsen (2013a) and refer to this metric as JP (note that Johannsen 2013a refers to this metric as the modified Kerr metric or MK).

In Boyer–Lindquist-like coordinates, the components for the JP metric are

$$\begin{eqnarray}\begin{array}{rcl}{g}_{{tt}}^{\mathrm{JP}} & = & -\displaystyle \frac{\tilde{{\rm{\Sigma }}}\left({\rm{\Delta }}-{a}^{2}{A}_{2}{\left(r\right)}^{2}{\sin }^{2}\theta \right)}{{\left[{A}_{1}(r)\left({a}^{2}+{r}^{2}\right)-{a}^{2}{A}_{2}(r){\sin }^{2}\theta \right]}^{2}},\\ {g}_{{rr}}^{\mathrm{JP}} & = & \displaystyle \frac{\tilde{{\rm{\Sigma }}}}{{A}_{5}(r){\rm{\Delta }}},\,\,\,{g}_{\theta \theta }^{\mathrm{JP}}=\tilde{{\rm{\Sigma }}},\\ {g}_{\phi \phi }^{\mathrm{JP}} & = & \displaystyle \frac{\tilde{{\rm{\Sigma }}}{\sin }^{2}\theta \left[{A}_{1}{\left(r\right)}^{2}{\left({a}^{2}+{r}^{2}\right)}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta \right]}{{\left[{A}_{1}(r)\left({a}^{2}+{r}^{2}\right)-{a}^{2}{A}_{2}(r){\sin }^{2}\theta \right]}^{2}},\\ {g}_{t\phi }^{\mathrm{JP}} & = & -\displaystyle \frac{a\tilde{{\rm{\Sigma }}}{\sin }^{2}\theta \left[{A}_{1}(r){A}_{2}(r)\left({a}^{2}+{r}^{2}\right)-{\rm{\Delta }}\right]}{{\left[{A}_{1}(r)\left({a}^{2}+{r}^{2}\right)-{a}^{2}{A}_{2}(r){\sin }^{2}\theta \right]}^{2}},\end{array}\end{eqnarray} \tag{ 7 }$$

where, to lowest order,

$$\begin{eqnarray}\begin{array}{rcl}{A}_{1}(r) & \equiv & \displaystyle \frac{{\alpha }_{13}{M}^{3}}{{r}^{3}}+1,\\ {A}_{2}(r) & \equiv & \displaystyle \frac{{\alpha }_{22}{M}^{2}}{{r}^{2}}+1,\\ {A}_{5}(r) & \equiv & \displaystyle \frac{{\alpha }_{52}{M}^{2}}{{r}^{2}}+1,\\ \tilde{{\rm{\Sigma }}} & \equiv & {\rm{\Sigma }}+\displaystyle \frac{{\epsilon }_{3}{M}^{3}}{r}.\end{array}\end{eqnarray} \tag{ 8 }$$

In this form, the metric has four free parameters, ₃, α₁₃, α₂₂, and α₅₂. The Kerr metric is recovered when all free parameters are set to zero. Johannsen (2013a) showed that the shape of the black hole shadow only depends on the parameters α₁₃ and α₂₂. Therefore, here, we will only vary those two parameters and set ₃ and α₅₂ to zero.

This metric describes a rotating vacuum spacetime that deviates from the Kerr metric but does not have pathologies outside of the event horizon for all spins below unity as long as the perturbation parameters satisfy (see also Figure 2 in Johannsen 2013a)

$$\begin{eqnarray}&&\begin{array}{l}{\alpha }_{13}\gt \displaystyle \frac{-{\left(M+\sqrt{{M}^{2}-{a}^{2}}\right)}^{3}}{{M}^{3}},\\ {\alpha }_{22}\gt \displaystyle \frac{-{\left(M+\sqrt{{M}^{2}-{a}^{2}}\right)}^{2}}{{M}^{2}}.\end{array}\end{eqnarray} \tag{ 9 }$$

For simplicity, in the present work we will only consider α₁₃ ≥ −1.0 and α₂₂ ≥ −1.0 since these are the lower limits that correspond to a maximally spinning black hole. Even though there is no upper limit on the values of these parameters, we simulate here shadows with α₁₃ ≤ 5 and α₂₂ ≤ 5. The event horizon for the JP metric does not depend on the perturbation parameters and coincides with the Kerr horizon,

$$\begin{eqnarray}&&{r}_{+}=M+\sqrt{{M}^{2}-{a}^{2}},\end{eqnarray} \tag{ 10 }$$

for all allowed values of the perturbation parameters.

2.3. The Modified Gravity Bumpy Kerr Metric

The metric proposed by Vigeland et al. (2011) uses bump functions to perturb the Kerr metric while still ensuring that the metric has three constants of motion. Here, we employ the metric in the form used by Gair & Yunes (2011) and Johannsen (2013b), which makes some simplifications to the original metric to ensure certain properties (see Gair & Yunes 2011 for details). We will employ the terminology in Johannsen (2013b) and refer to this metric as the modified gravity bumpy Kerr metric or MGBK. The MGBK metric is defined perturbatively in covariant form such that

$$\begin{eqnarray}&&{g}_{\mu \nu }^{\mathrm{MGBK}}={g}_{\mu \nu }^{{\rm{K}}}+{h}_{\mu \nu }^{\mathrm{MGBK}},\end{eqnarray} \tag{ 11 }$$

where the perturbations to each metric component are

$$\begin{eqnarray}&&\begin{array}{l}{h}_{{tt}}^{\mathrm{MGBK}}=-\displaystyle \frac{{{aP}}_{2}}{{{MP}}_{1}}{h}_{t\phi }-\,\displaystyle \frac{a{\rm{\Delta }}{{\rm{\Sigma }}}^{2}}{2{{MP}}_{1}}\displaystyle \frac{\partial {h}_{t\phi }}{\partial r}-\displaystyle \frac{{a}^{2}{{\rm{\Delta }}}^{2}{\rm{\Sigma }}{\sin }^{2}\theta }{2{{MP}}_{1}}\displaystyle \frac{d{\bar{\gamma }}_{1}}{{dr}}-\,\displaystyle \frac{a{{\rm{\Delta }}}^{2}{\sin }^{2}\theta \left(2{a}^{2}{Mr}{\sin }^{2}\theta +\hat{{\rm{\Sigma }}}\right)}{2{{MP}}_{1}}\displaystyle \frac{d{\bar{\gamma }}_{3}}{{dr}}-\,\displaystyle \frac{{a}^{2}{{\rm{\Delta }}}^{2}{\sin }^{2}\theta ({\rm{\Sigma }}-4{Mr})}{2{{MP}}_{1}}\displaystyle \frac{d{\bar{\gamma }}_{4}}{{dr}}\\ \,+\,\displaystyle \frac{2{a}^{2}{\rm{\Delta }}{\bar{\gamma }}_{1}{r}^{2}{\sin }^{2}\theta }{{P}_{1}}+\displaystyle \frac{{\rm{\Delta }}{\bar{\gamma }}_{1}{\hat{\rho }}^{2}\left({a}^{2}+{r}^{2}\right)}{{P}_{1}}-\,\displaystyle \frac{a{\rm{\Delta }}{\bar{\gamma }}_{3}{P}_{3}{\sin }^{2}\theta }{M{\rm{\Sigma }}{P}_{1}}+\displaystyle \frac{2{\rm{\Delta }}{\bar{\gamma }}_{4}{P}_{4}}{{\rm{\Sigma }}{P}_{1}},\\ {h}_{{rr}}^{\mathrm{MGBK}}=-\displaystyle \frac{{\rm{\Sigma }}{\bar{\gamma }}_{1}}{{\rm{\Delta }}},\\ {h}_{t\phi }^{\mathrm{MGBK}}=\,-\displaystyle \frac{{\sin }^{2}\theta \left[{a}^{2}M{\gamma }_{\mathrm{3,1}}+{{aM}}^{2}({\gamma }_{\mathrm{1,2}}+{\gamma }_{\mathrm{4,2}})+{M}^{3}{\gamma }_{\mathrm{3,3}}\right]}{{\rm{\Sigma }}},\\ {h}_{\phi \phi }^{\mathrm{MGBK}}=-\displaystyle \frac{{\left({a}^{2}+{r}^{2}\right)}^{2}}{{a}^{2}}{h}_{{tt}}+\displaystyle \frac{{\rm{\Delta }}{\rm{\Sigma }}{\bar{\gamma }}_{1}}{{a}^{2}}-\,\displaystyle \frac{2\left({a}^{2}+{r}^{2}\right)}{a}{h}_{t\phi }-\displaystyle \frac{2{{\rm{\Delta }}}^{2}{\bar{\gamma }}_{3}{\sin }^{2}\theta }{a}+\displaystyle \frac{2{{\rm{\Delta }}}^{2}{\bar{\gamma }}_{4}}{{a}^{2}},\end{array}\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{P}_{1} & = & -{a}^{4}\left({r}^{2}-{a}^{2}\right){\cos }^{4}\theta -2{a}^{2}\left({a}^{4}+{a}^{2}{r}^{2}+2{r}^{4}\right)\;{\cos }^{2}\theta +{r}^{2}\left(2{a}^{4}+5{a}^{2}{r}^{2}+{r}^{4}\right),\\ {P}_{2} & = & -{a}^{4}{\cos }^{4}\theta (r-M)-{a}^{2}{\cos }^{2}\theta \,\left[{r}^{2}(2M+r)-{a}^{2}(r-2M)\right]+\,{r}^{2}\left(2{a}^{2}M+{a}^{2}r+{{Mr}}^{2}\right),\\ {P}_{3} & = & {a}^{4}{\rm{\Delta }}{\cos }^{4}\theta (r-2M)+2{a}^{2}{r}^{3}{\cos }^{2}\theta \,\left({a}^{2}+4{M}^{2}-2{Mr}+{r}^{2}\right)+{r}^{3}(r-2M)\,\left[{a}^{2}(4M+r)+{r}^{2}(2M+r)\right],\\ {P}_{4} & = & -2{a}^{2}{r}^{2}{\cos }^{2}\theta \left(-{a}^{2}-2{Mr}+{r}^{2}\right)+{r}^{3}\left[{a}^{2}(3r-4M)+{r}^{2}(r-2M)\right]-{a}^{4}{\cos }^{4}\theta \left({a}^{2}-2{Mr}+3{r}^{2}\right),\end{array}\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{\hat{\rho }}^{2} & = & {r}^{2}-{a}^{2}{\cos }^{2}\theta ,\\ \hat{{\rm{\Sigma }}} & = & {\left({a}^{2}+{r}^{2}\right)}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta ,\\ {\bar{\gamma }}_{1} & = & \displaystyle \frac{{\gamma }_{\mathrm{1,2}}{M}^{2}}{{r}^{2}},\\ {\bar{\gamma }}_{3} & = & \displaystyle \frac{{\gamma }_{\mathrm{3,1}}{{Mr}}^{2}+{\gamma }_{\mathrm{3,3}}{M}^{3}}{{r}^{4}},\\ {\bar{\gamma }}_{4} & = & \displaystyle \frac{{\gamma }_{\mathrm{4,2}}{M}^{2}}{{r}^{2}}.\end{array}\end{eqnarray} \tag{ 14 }$$

In the equations above, Σ and Δ are the same as in the Kerr metric (see Equation (2)). Here we have resolved a typographical error present in the definition of P₁ in Vigeland et al. (2011), which was also corrected in Vigeland (2012) but not utilized in subsequent literature. The resolution of this error allows us to explore a wider range of values for the free parameters of this metric, γ_3,1, γ_3,3, γ_1,3, and γ_4,2, than was done in Johannsen (2013b). We allow all free parameters to vary between 0.0 and 2.0. In this form, the event horizon of this metric coincides with the Kerr horizon. This metric is also not Ricci flat.

We simulate a large set of black hole shadows that probe the allowed parameter space of the metrics described in the previous section. We perform ray-tracing simulations using the algorithm of Psaltis & Johannsen (2012), where we solve the geodesic equations in each metric to derive the black hole shadows.

All of the metrics are axisymmetric, which allows us to use the two associated killing vectors to simplify the geodesic equations. We solve the simplified geodesic equations with a fourth-order Runge–Kutta integration scheme, with adaptive step size. We integrate the photon trajectories backwards starting at the observer's image plane and ending when the photon either comes within 10% of the event horizon, or escapes to infinity. We denote by $x^{\prime}$ and $y^{\prime}$ the coordinates on the image plane (see Figure 1 in Johannsen & Psaltis 2010b for the geometry used for these simulations). For the QK metric we exclude the region inside $r=2.48\,M$ to avoid encountering any pathologies.

We describe the shape of a black hole shadow in terms of its image-plane radius $R\equiv \sqrt{x{{\prime} }^{2}+y{{\prime} }^{2}}$ at different orientation angles $\psi ={\cos }^{-1}(x^{\prime} /R)$. We place 200 evenly spaced resolution elements along ψ, for 0 ≤ ψ ≤ 2π, and, for each value of this orientation angle, use the bisection method to find the boundary between photon trajectories that trace back into the black hole horizon and those that escape to infinity, i.e., the black hole shadow; we terminate the bisection at an interval ≲10⁻⁶ M. Since we use Boyer–Lindquist-like coordinates for all metrics, there is a pole that coincides with the spin axis of the black hole that often introduces numerical errors in the calculations. For this reason, we exclude orientations that are within 0.001 rad of ψ = π/2 or ψ = 3π/2 and interpolate between the nearby data points within these intervals.

In Figure 2, we show how the shadow that results from the QK metric depends on the quadrupole deviation parameter _QK, while keeping the spin and the inclination angle of the observer constant. As discussed in Johannsen & Psaltis (2010b), the quadrupole deviation parameter in this metric introduces asymmetry to the shape of the shadow. However, because we opted here not to expand the geodesic equations in this metric to first order in _QK, we can only consider small values for this parameter before encountering metric pathologies. This results in shadow shapes that show only small deviations from the Kerr shadows.

Zoom In Zoom Out Reset image size

In Figure 3, we show how the shadows that result from the JP metric depend on the deviation parameters α₁₃ and α₂₂. The effect of changing the parameter α₂₂ on the shadow is similar to increasing the spin. For a Kerr metric, the shadow is approximately circular up to a_* ≈ 0.9, after which a perturbation on the left side of the shadow appears (see Figure 1). For the JP metric, increasing α₂₂ increases the perturbation to the left side of the shadow far beyond the maximum for a Kerr black hole.

Zoom In Zoom Out Reset image size

On the other hand, the effect of changing the parameter α₁₃ on the shadow is similar to increasing the black hole mass. Kerr black hole shadows have radii between 4.8M and 5.2M for all spins. For the JP metric, increasing the parameter α₁₃ tends to increase the size of the shadow, while making it more circular. Within the set of shadows we explored, the radius range for the JP metric is 5M ≲ R ≲ 6M, indicating that a JP black hole with α₁₃ = 5 can have a black hole shadow that is ≈20% larger than a Kerr black hole of equivalent mass.

In Figure 4 we show how the shadows that result from the MGBK metric depend on the perturbation parameters γ_3,1, γ_3,3, γ_1,2, and γ_4,2. The perturbations created by the four free parameters of the MGBK metric have some similarities. However, unlike the case of the other metrics, the shadows can look qualitatively different at low versus high a_*. For all parameters, the shadows seem to have "dimples" along the $x^{\prime} =0$ line.

Zoom In Zoom Out Reset image size

When a_* is small, as is shown in the top right panel of Figure 4, the shadows are relatively left–right symmetric and can become significantly elongated in the horizontal direction. The parameter γ_3,3 creates the most extreme shadows of this kind, as shown in the figure. In contrast, γ_3,1 has a negligible effect when a_* is small.

When a_* is large, as is shown in the top left and bottom panels, the shadows become asymmetric with the right side becoming larger than the left due to the effect of frame dragging, as in the Kerr shadows in Figure 1. The parameters γ_3,1 and γ_3,3 do not affect the overall size of the shadow for large a_*, while the parameters γ_1,2 and γ_4,2 can create larger shadows for large a_* as shown. For high a_* and high γ_1,2 shadows can become flat on the left side, similar to the behavior of Kerr shadows for very high spin.

Frame-dragging effects cause a significant fraction of the shadows we have calculated to be displaced from the image origin that is centered on the black hole. Moreover, the effect of some perturbation parameters is simply to rescale the shadow sizes. Because our goal is to treat the set of shadows as an ensemble and compactly represent their shapes, we will first remove all trivial transformations. A simple translation or rescaling can be easily accounted for after the fact and does not need to be included as part of the PCA. By rescaling and recentering the shadows before analyzing their shapes, we can significantly reduce the number of parameters that will later be needed to fit to data. We emphasize that we do not wish to discard the information about the original size and displacement of the shadows but rather choose to consider it separately.

There are many possible ways of defining the center of a perturbed circle, each resulting in a different set of parameterizations for the shadow shapes. For this work, we experimented with various methods and chose the method that resulted in the least number of PCA components (see Section 5) necessary for accurate reconstructions. We also preferred a procedure for standardizing shadows that could be easily applied to the data as well as the simulations. We chose to recenter and rescale each shadow such that it is as large as possible while still being contained within a circle centered at the origin with a radius equal to $\sqrt{27}M$, i.e., the radius of the shadow of a Schwarzschild black hole.

We denote the amount by which each shadow is displaced in the $x^{\prime}$ direction, i.e., perpendicular to the spin axis, by x₀, such that a shadow that is displaced to the right will have a positive x₀. Almost all of our shadows have negligible displacement in the $y^{\prime}$ direction, i.e., along the spin axis, so we do not include a vertical displacement in this discussion.

Before recording the shape of each shadow in terms of the image-plane radius R along different orientation angles ψ (see Section 3), we first recenter it by subtracting (x₀, y₀) from all coordinate pairs. We then define the scale of the shadow R_scale as the largest radius along the perimeter of the shadow and then rescale the shadows such that the largest radius becomes equal to $\sqrt{27}M$.

We further define the asymmetry (A_asym) of a shadow as

$$\begin{eqnarray}&&{A}_{\mathrm{asym}}=\displaystyle \frac{{\sigma }_{R}}{\bar{R}},\end{eqnarray} \tag{ 15 }$$

where σ_R is the standard deviation and $\bar{R}$ is the mean of the radius of each shadow along different orientation angles ψ.

In Figures 5–8 we show the dependence of the displacement (x₀), the scaling factor (R_scale), and the asymmetry (A_asym) of the shadows on the perturbation parameters of the various metrics and the black hole spin; for comparison, we also show the same parameters for the Kerr metric. As expected, for all metrics, the shadows are displaced to the right with increasing black hole spin, with the perturbation parameters having only a secondary effect. At the same time, changing most of the perturbation parameters introduces changes in the overall scale R_scale of the shadows to amounts larger than the ±4% range obtained for the Kerr metric. Both the JP and the MGBK metrics can create asymmetries, A_asym, that are much larger than what is seen for the Kerr metric.

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

Before performing PCA on the recentered and rescaled set of shadows, we interpolate each shadow so that they are all evaluated on evenly spaced values of ψ and aligned with all other shadows. In Figure 9 we show the full set of shadows after they have been rescaled and recentered. A total of 23,887 shadows are included in the set: 15,625 for the MGBK metric, 7605 for the JP metric, 522 for the QK metric, and 135 for the Kerr metric. For the QK metric, we allowed the spin of the black hole to vary between 0.0 and 0.4 in intervals of 0.05 and the perturbation parameter _QK to vary between 0.0 and 0.35 in intervals of 0.05. (As previously noted, the limit of allowed values of the perturbation parameter _QK is smaller for certain combinations of a_* and i.) For the JP metric, we allowed the spin of the black hole to vary between 0.1 and 0.9 in intervals of 0.2, and both of the perturbation parameters, α₁₃ and α₂₂, to vary between −1 and 5 in intervals of 0.5. For the MGBK metric we allowed a_* to vary from 0.1 to 0.9 in intervals of 0.2 and allowed the perturbation parameters, γ_3,1, γ_3,3, γ_1,2, and γ_4,2, to vary between 0.0 and 2.0 in intervals of 0.5. For the Kerr metric we allowed a_* to vary between 0.1 and 0.9 in intervals of 0.1 and also included the values 0.95, 0.99, 0.9999, 0.99999, 0.999999, and 0.99999999. For all metrics except MGBK, we allowed the inclination angle i to vary between 10° and 90° in intervals of 10°; for MGBK we used intervals of 20°.

Zoom In Zoom Out Reset image size

In order to compare the shadow shapes we have simulated to observations, we would like to express their shapes using a small number of parameters that can be directly inferred from the data. Abdujabbarov et al. (2015) followed a similar approach by representing the shadows as a sum of Legendre polynomials. Here we explore the use of PCA to compactly represent the space spanned by Kerr and non-Kerr shadows.

The PCA algorithm diagonalizes the covariance matrix of the data to find a set of orthogonal basis vectors (principal components) that are ordered such that the first principal component accounts for the largest possible percentage of the variance in the data set, the second principal component accounts for the second largest percentage of the variance, and so on (see, e.g., Boroson & Lauer 2010). For data sets that are correlated, this can result in significant dimensionality reduction since it may be possible to reconstruct the original data set using only a handful of PCA components.

We represent each shadow as a column vector S_n ≡ S_in, which corresponds to the radius of the ith orientation point of the nth shadow. In addition to rescaling and recentering as was discussed above, we also subtract from each shadow shape the constant $\sqrt{27}M$, such that all circular shadows would have a radius of zero after rescaling, recentering, and subtraction.

We define the covariance matrix C as

$$\begin{eqnarray}&&C=\displaystyle \frac{1}{m}\sum _{n=1}^{m}({{\boldsymbol{S}}}_{n}-\sqrt{27}M){\left({{\boldsymbol{S}}}_{n}-\sqrt{27}M\right)}^{T},\end{eqnarray} \tag{ 16 }$$

and find a basis of eigenvectors, or "eigenshadows," by diagonalizing C such that

$$\begin{eqnarray}&&C{{\boldsymbol{u}}}_{k}={\lambda }_{k}{{\boldsymbol{u}}}_{k},\end{eqnarray} \tag{ 17 }$$

where u_k is the kth eigenshadow, λ_k is the kth eigenvalue, and k = 1,2, ... , m. Although PCA returns m eigenshadows, they are ordered such that only the first few eigenshadows are needed to reconstruct the majority of the variance within the data set. The associated eigenvalue of an eigenshadow indicates the percentage of the variance within the data set that this particular eigenshadow accounts for. Because of this, we normalize each eigenshadow such that

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{k}^{2}={\lambda }_{k},\end{eqnarray} \tag{ 18 }$$

and normalize the eigenvalues such that

$$\begin{eqnarray}&&\sum _{k=1}^{m}{\lambda }_{k}=1.\end{eqnarray} \tag{ 19 }$$

Once we have derived a basis of eigenshadows, we can reconstruct each shadow as a linear combination of the eigenshadows. Specifically, we write

$$\begin{eqnarray}&&{{\boldsymbol{S}}}_{n}=\sqrt{27}M+\sum _{k=1}^{m}{a}_{{nk}}{{\boldsymbol{u}}}_{k},\end{eqnarray} \tag{ 20 }$$

where the amplitude a_nk is the coefficient multiplying the kth eigenshadow in the linear combination of the nth shadow. The overall sign of each eigenshadow is arbitrary and any eigenshadow can have both positive and negative coefficients (a_nk). For simplicity, we choose the sign of each eigenshadow such that the amplitude with the largest magnitude is positive. In the next section we will also discuss the relative contribution of eigenshadow k to a shadow n, which we define as

$$\begin{eqnarray}&&{a}_{{nk}}^{{\prime} }=\displaystyle \frac{{a}_{{nk}}{\sigma }_{{{\boldsymbol{u}}}_{k}}}{\sqrt{27}M},\end{eqnarray} \tag{ 21 }$$

where ${\sigma }_{{{\boldsymbol{u}}}_{k}}$ is the standard deviation of the kth eigenshadow defined as

$$\begin{eqnarray}&&{\sigma }_{{{\boldsymbol{u}}}_{k}}=\sqrt{\displaystyle \frac{1}{N}\sum _{i=1}^{N}{\left({u}_{{ik}}-\langle {{\boldsymbol{u}}}_{k}\rangle \right)}^{2}}.\end{eqnarray} \tag{ 22 }$$

Here, u_ik is the ith element of the u_k eigenshadow, and $\langle {{\boldsymbol{u}}}_{k}\rangle$ is the average of the elements of the u_k eigenshadow. We will also refer to ${a}_{{nk}}^{{\prime} }$ as the rescaled amplitude.

In principle, to reconstruct a shadow to high accuracy, we would need all m eigenshadows. However, because PCA components are ordered by their relative importance, the first few eigenshadows can already reconstruct the large scale features of each shadow, as long as the eigenvalues decline reasonably fast. This allows for significant dimensionality reduction as only a handful of eigenshadows may be needed depending on the desired accuracy. In the following sections, we will explore how many eigenshadows are needed for accurate reconstructions.

In Figure 10, we show the full set of shadows unwrapped into curves of radius R as a function of the orientation angle ψ. We also show the first six PCA components that were derived from the set of shadows. In this figure, we have multiplied each PCA component by the largest amplitude needed to reconstruct the shadows within the set, effectively showing the largest possible contribution that each eigenshadow has to any shadow within the set.

Zoom In Zoom Out Reset image size

In Figure 11, we show the spectrum of PCA eigenvalues. The rapid decrease in PCA eigenvalues indicates that only a small number of PCA components are needed to reproduce the majority of the variance in the data set. Specifically, the first PCA component accounts for 81.6% of the variance in the data set, the second component accounts for 13.9%, and only four PCA components are needed to account for 98.6% of the variance.

Zoom In Zoom Out Reset image size

In this section we quantify the accuracy of the PCA reconstructions using only a small number of components and contrast it to reconstructions using Legendre polynomials (as done in, e.g., Abdujabbarov et al. 2015). Additionally, we explore the relation between the amplitudes of the PCA components needed to reconstruct a particular shadow and the parameters of the underlying metrics.

We define the maximum fractional difference (ξ) for the PCA reconstruction of each shadow as

$$\begin{eqnarray}&&\xi =\displaystyle \frac{\max (| {S}_{\mathrm{in}}-{S}_{\mathrm{in}}^{{\prime} }| )}{\sqrt{27}M},\end{eqnarray} \tag{ 23 }$$

where S_in denotes the complete shadow and ${S}_{\mathrm{in}}^{{\prime} }$ denotes the reconstructed shadow. We apply this metric to reconstructions with either PCA components or with Legendre polynomials. For consistency with the PCA reconstructions we perform the Legendre polynomial reconstructions on the recentered and rescaled shadows.

In Figure 12 we plot the maximum fractional difference as a function of the asymmetry of each shadow for reconstructions using 3–6 PCA components and 3–6 Legendre polynomials. As expected, both methods have similar performance for shadows with small asymmetry. However, the maximum fractional difference in the Legendre polynomial reconstructions is significantly higher than that of the PCA reconstruction for shadows with higher asymmetry. A reconstruction with only five Legendre polynomials for shadows with a ∼20% asymmetry will lead to a fractional difference as high as 25%. In contrast, only five PCA components are needed to reconstruct all shadows within the ensemble with ≲5% maximum fractional error. Adding a sixth PCA component does not substantially improve the reconstruction.

Zoom In Zoom Out Reset image size

The shadows that have the highest (≳5%) fractional maximum error in reconstructions are MGBK shadows with high perturbation parameters and high values of spin, such as those shown in the bottom panels of Figure 4. Most of these shadows have small "dimples" at $x^{\prime} =0$ before they are aligned, which are displaced to different x locations, depending on spin, after alignment. These dimples make the shadows particularly hard to reconstruct because their varying locations cannot be captured by a small number of PCA components.

In Figure 13 we show two sample reconstructions for shadows with high asymmetry using either four PCA components or four Legendre polynomials. The left panel corresponds to the MGBK metric with parameters a_* = 0.1, i = 90°, γ_3,1 = 0.0, γ_3,3 = 2.0, γ_1,2 = 0.0, and γ_4,2 = 2.0 and the right panel corresponds to the JP metric with parameters a_* = 0.9, i = 70°, α₁₃ = 2.5, and α₂₂ = 5.0. Only a small number of PCA components are needed to generate high fidelity reconstructions. Reconstructions with the same number of Legendre polynomials result in significantly higher errors.

Zoom In Zoom Out Reset image size

We have calculated an ensemble of over 20,000 simulated black hole shadows that probe the allowed space of several parameterized non-Kerr metrics. We applied PCA to our ensemble of shadows and derived an orthogonal basis for shadow shapes. We then used this analysis to reduce the dimensionality of the set of eigenshadows and showed that only five PCA components are required to reproduce all shadows within the ensemble, with a maximum fractional difference (ξ) of ≲5%. We compared the accuracy of the PCA reconstructions to reconstructions using Legendre polynomials and found that, for shadows with high asymmetry, the PCA reconstructions introduce errors that are nearly 10 times lower than the Legendre polynomial reconstructions.

This approach allows us to measure or place constraints on deviations from the Kerr metric using black hole images. Given a set of measurements and choosing a particular metric, we can convert the measured coefficients of the PCA components to metric parameters. Figures 14–17 show sample mappings between the rescaled amplitudes (${a}_{{nk}}^{{\prime} }$) of the PCA reconstructions of shadows and the parameters of the underlying metrics. Such mapping will allow us to place metric-specific constraints on deviations from the Kerr metric.

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

The accuracy with which metric-specific constraints can be imposed and the ability to break degeneracies between different deviation parameters will depend on both the availability of prior information on the properties of the observed black hole and on the interferometric coverage of the observations. For example, as Figure 3 shows, in the JP metric, a nearly circular black hole shadow can be obtained either for small spin and minimal deviations from Kerr or large spin and a large value for the deviation parameter α₁₃. The two cases, however, result in substantially different shadow sizes, which will correspond to different underlying black hole masses. Prior knowledge of the mass of the black hole, e.g., from monitoring the orbits of stars around it, as is the case for Sgr A*, will help in breaking such a degeneracy. Similarly, nearly circular shadows can also be obtained by reducing the inclination of the observer with respect to the spin axis of the black hole. Prior information of this inclination, as is the case for M87, will be valuable in breaking such degeneracies.

The current interferometric coverage available for EHT observations allows for measuring only global properties of the black hole shadows, such as their sizes and overall asymmetries, and not for tracing out their precise shapes (Event Horizon Telescope Collaboration et al. 2019b). As Figures 7 and 8 show, such general measurements can be used primarily to place correlated constraints on deviation parameters for different metrics and not to measure these parameters independently. The addition of more EHT baselines on Earth and, in the future, a more complete interferometric coverage using rapidly orbiting, space-based stations will allow measurement of the precise shapes of black hole shadows and reduce possible correlations between parameters that describe deviations from the Kerr metric.

A set of measured PCA amplitudes of the shadow shape may also allow us to distinguish between different metrics. Figure 18 shows the distribution of amplitudes of the first two PCA components for each metric we considered here. Measuring a large positive ratio between the second and the first PCA amplitudes will point toward the Kerr metric or modifications similar to those described within the JP metric. On the other hand, measuring a negative ratio between the second and the first PCA amplitudes will point toward modifications similar to those described in the MGBK metric.

Zoom In Zoom Out Reset image size

The eigenshadows we derived can be directly applied to the results of EHT observations (see, e.g., Event Horizon Telescope Collaboration et al. 2019b) as they can be incorporated into a Hough/Radon transform and utilized in an edge detection algorithm (see Psaltis et al. 2015). The location of the shadow in the image plane and the size of the shadow (R_scale) can also be incorporated into the Hough/Radon transform and compared to the R_scale and x₀ of the shadows in our data set. The outcome of this application will be a measurement of (or a constraint on) the coefficients of the various PCA components. This measurement will be agnostic of the metric of the compact object.

We thank C.K. Chan, T. Lauer, Z. Younsi, C. Raithel, D. Ball, J. Kim, T. Trent, M. Rose, and M. de Laurentis for useful discussions. L.M. gratefully acknowledges support from NSF GRFP grant DGE 1144085, NSF PIRE grant 1743747, and an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award No. AST-1903847. D.P. and F.O. gratefully acknowledge support from NSF PIRE grant 1743747, NSF AST-1715061, and Chandra Award No. TM8-19008X for this work. All ray tracing and PCA calculations were performed with the El Gato GPU cluster at the University of Arizona that is funded by NSF award 1228509.