Testing General Relativity with X-ray reflection spectroscopy: The Konoplya-Rezzolla-Zhidenko parametrization

TL;DR

This paper explores the use of X-ray reflection spectroscopy to test deviations from the Kerr black hole model in general relativity, implementing a parametrized non-Kerr metric and deriving observational constraints.

Contribution

It introduces the implementation of the Konoplya-Rezzolla-Zhidenko black hole metric in the extsc{relxill extunderscore nk} framework for testing non-Kerr black holes with X-ray data.

Findings

First X-ray constraints on deformation parameters of the non-Kerr metric.

Analysis of how strong-field parameters affect observable spectra.

Comparison of X-ray constraints with gravitational wave and black hole imaging results.

Abstract

X-ray reflection spectroscopy is a promising technique for testing general relativity in the strong field regime, as it can be used to test the Kerr black hole hypothesis. In this context, the parametrically deformed black hole metrics proposed by Konoplya, Rezzolla \& Zhidenko (Phys. Rev. D93, 064015, 2016) form an important class of non-Kerr black holes. We implement this class of black hole metrics in \textsc{relxill\_nk}, which is a framework we have developed for testing for non-Kerr black holes using X-ray reflection spectroscopy. We perform a qualitative analysis of the effect of the leading order strong-field deformation parameters on typical observables like the innermost stable circular orbits and the reflection spectra. We also present the first X-ray constraints on some of the deformation parameters of this metric, using \textit{Suzaku} data from the supermassive black hole in Ark~564, and compare them with those obtained (or expected) from other observational techniques like gravitational waves and black hole imaging.

Tables

Table 1. Table 1: Summary of the best-fit values for the spectral models 1 and 2, the latter for each δ i subscript 𝛿 𝑖 \delta_{i} case. Uncertainties are reported at 90% confidence level. All values are reported up to two significant digits after the decimal, except a ∗ subscript 𝑎 a_{*} and δ i subscript 𝛿 𝑖 \delta_{i} which are reported to three significant digits after the decimal. P P \mathrm{P} indicates the relevant parameter being unbounded up to its explored range and ϵ italic-ϵ \epsilon indicates uncertainty too small up to the precision reported here. See the text for more details.

Model	1	2$\left({\delta }_1\right)$	2$\left({\delta }_2\right)$	2$\left({\delta }_3\right)$	2$\left({\delta }_4\right)$	2$\left({\delta }_5\right)$	2$\left({\delta }_6\right)$
tbabs
$N_H/{10}^{20}c{m}^{-2}$	${6.74}^{\ast }$	${6.74}^{\ast }$	${6.74}^{\ast }$	${6.74}^{\ast }$	${6.74}^{\ast }$	${6.74}^{\ast }$	${6.74}^{\ast }$
zpowerlaw
$\Gamma$	${2.78}_{-0.01}^{+ϵ}$	${2.86}_{-0.04}^{+0.05}$	${2.86}_{-0.04}^{+0.06}$	${2.86}_{-0.04}^{+0.06}$	${2.85}_{-0.03}^{+0.01}$	${2.87}_{-0.03}^{+0.04}$	${2.87}_{-0.02}^{+0.05}$
$N$	$0.68$	${0.39}_{-0.04}^{+0.03}$	${0.39}_{-0.04}^{+0.02}$	${0.39}_{-0.06}^{+0.03}$	${0.39}_{-0.03}^{+0.02}$	${0.39}_{-0.03}^{+0.02}$	${0.39}_{-0.04}^{+0.02}$
relxill_nk
$q$	-	${8.20}_{-0.05}^{+P}$	${8.66}_{-0.05}^{+P}$	${8.57}_{-0.49}^{+0.79}$	${8.52}_{-0.36}^{+0.67}$	${8.83}_{-0.63}^{+0.44}$	${8.39}_{-0.36}^{+0.48}$
$a_{\ast }$	-	${0.998}_{-0.022}^{+P}$	${0.998}_{-0.010}^{+P}$	${0.998}_{-0.012}^{+P}$	${0.986}_{-0.022}^{+0.009}$	${0.986}_{-0.003}^{+0.006}$	${0.986}_{-0.013}^{+0.005}$
$i$ [deg]	-	${27}_{-P}^{+2}$	${27}_{-P}^{+3}$	${27}_{-P}^{+P}$	${27}_{-P}^{+4}$	${29}_{-P}^{+P}$	${27}_{-P}^{+P}$
$\log \xi$	-	${3.00}_{-0.53}^{+0.12}$	${3.00}_{-0.18}^{+0.07}$	${3.00}_{-0.24}^{+0.09}$	${3.00}_{-0.16}^{+0.04}$	${3.00}_{-0.16}^{+0.03}$	${2.98}_{-0.11}^{+0.08}$
$A_{\text{Fe}}$	-	${0.83}_{-0.18}^{+0.13}$	${0.84}_{-0.09}^{+0.14}$	${0.85}_{-0.18}^{+0.14}$	${0.85}_{-0.15}^{+0.07}$	${0.81}_{-0.09}^{+0.09}$	${0.88}_{-0.09}^{+0.09}$
$\delta$-type	-	${\delta }_1$	${\delta }_2$	${\delta }_3$	${\delta }_4$	${\delta }_5$	${\delta }_6$
$\delta$-value	-	${0.197}_{-0.464}^{+0.087}$	$-{0.254}_{-0.117}^{+0.469}$	$-{0.414}_{-0.035}^{+1.452}$	${2.928}_{-3.475}^{+P}$	$-{0.896}_{-P}^{+0.818}$	$-{0.610}_{-4.262}^{+0.510}$
$N/{10}^{-2}$	-	${3.56}_{-0.52}^{+3.61}$	${3.51}_{-0.28}^{+0.84}$	${3.51}_{-0.44}^{+3.06}$	${3.28}_{-0.33}^{+0.50}$	${3.91}_{-0.28}^{+0.31}$	${3.85}_{-0.28}^{+0.31}$
xillver
$N/{10}^{-3}$	-	${5.20}_{-1.64}^{+1.32}$	${5.26}_{-1.59}^{+1.80}$	${5.25}_{-1.66}^{+1.83}$	${5.28}_{-1.74}^{+1.59}$	${5.32}_{-1.80}^{+1.62}$	${5.16}_{-1.80}^{+1.62}$
${\chi }^2/dof$	$3.30$	$1.10$	$1.10$	$1.10$	$1.10$	$1.10$	$1.10$

Table 2. Table 2: Lower and upper bounds on the KRZ deformation parameter δ 1 subscript 𝛿 1 \delta_{1} from three different techniques. In the case of gravitational waves, the best constraints are reported, as obtained from the event GW191226. In the other two cases, only one source/observation has been used as of now to estimate the constraint.

Lower bound	Upper bound	Source	Technique
–0.16	+0.17	GW151226	Gravitational wave inspiral Cardenas-Avendano et al. (2020a)
–0.27	+0.28	Ark 564	This work
–1.2	+1.3	M${87}^{\ast }$	BH shadow size Psaltis et al. (2020)

Equations

$$d{s}^2$$

$$B=1+\sum _{i=0}^{\infty }\left\{b_{i0}(1-x)+{\tilde{B}}_i(x)(1-x^2)\right\}{y}^i,$$

$${\tilde{B}}_i(x)=\frac{b_{i1}}{1+\frac{b_{i2}x}{1+\frac{b_{i3}x}{1+\cdots }}}.$$

$$N^2$$

$$B$$

$$W$$

$$K^2$$

$$\begin{array}{r}{r}_0=1+\sqrt{1-a_{\ast }^2},{ϵ}_0=\frac{2-r_0}{r_0},a_{20}=\frac{2a_{\ast }^2}{r_0^3},a_{21}=-\frac{a_{\ast }^4}{r_0^4}+{\delta }_6,\\ k_{00}=k_{22}=k_{23}=\frac{a_{\ast }^2}{r_0^2},k_{20}=0,k_{21}=\frac{a_{\ast }^4}{r_0^4}-\frac{2a_{\ast }^2}{r_0^3}-{\delta }_6,w_{00}=\frac{2a_{\ast }}{r_0^2},\end{array}$$

$${\delta }_1$$

$${\delta }_1>\frac{4r_0-3r_0^2-a^2}{r_0^2},$$

$$\displaystyle\delta_{2},\delta_{3}\left\{\begin{array}[]{l}>\\ <\\ \end{array}-\frac{4}{a^{3}}(1-\sqrt{1-a^{2}})\quad\begin{array}[]{l}\rm{if}\;a>0\\ \rm{if}\;a<0,\end{array}\right.$$

$${\delta }_4,{\delta }_5>-1,$$

$${\delta }_6<\frac{r_0^2}{4-a^2}.$$

$$F_o({\nu }_o)=\int I_o({\nu }_o,X,Y)d\tilde{\Omega }=\frac{1}{D^2}{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}{\int }_0^1\pi r_e\frac{g^2}{\sqrt{g^{\ast }(1-g^{\ast })}}f(g^{\ast },r_e,i)I_e({\nu }_e,r_e,{\vartheta }_e)d{g}^{\ast }d{r}_e.$$

$$f(g^{\ast },r_e,i)=\frac{1}{\pi r_e}g\sqrt{g^{\ast }(1-g^{\ast })}\left|\frac{\partial \left(X,Y\right)}{\partial \left(g^{\ast },r_e\right)}\right|.$$

$$g^{\ast }=\frac{g-g_{\min }}{g_{\max }-g_{\min }},$$

$$F_o({\nu }_o)=\frac{1}{D^2}{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}{\int }_0^1\pi r_e\frac{g^2}{\sqrt{g^{\ast }(1-g^{\ast })}}{f}^{(1)}(g^{\ast },r_e,i)I_e({\nu }_e,r_e,{\vartheta }_e^{(1)})d{g}^{\ast }d{r}_e$$

$$+\frac{1}{D^2}{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}{\int }_0^1\pi r_e\frac{g^2}{\sqrt{g^{\ast }(1-g^{\ast })}}{f}^{(2)}(g^{\ast },r_e,i)I_e({\nu }_e,r_e,{\vartheta }_e^{(2)})d{g}^{\ast }d{r}_e.$$

$$\begin{array}{rl}-0.27<& {\delta }_1<0.28,\\ -0.37<& {\delta }_2<0.22,\\ -0.45<& {\delta }_3<1.04,\\ -0.55<& {\delta }_4,\\ -1<& {\delta }_5<-0.08(0.45),\\ -4.9<& {\delta }_6<-0.10(0.36).\end{array}$$

Full text

I Introduction

Einstein’s theory of gravity, since its proposition over a century ago, has been applied to a variety of astrophysical phenomena in our Universe. Over these years, it has emerged as the standard framework for describing the spacetime in the presence of gravitational objects. While largely successful in the weak-field tests Will (2014), only recently the strong-field predictions of Einstein’s gravity (the general theory of relativity, GR hereafter) have become testable in a variety of ways Bambi (2017a); Abbott et al. (2019); Akiyama et al. (2019a). Presence of a zoo of alternative theories of gravity, which address shortcomings of GR with respect to observations e.g., dark matter and dark energy, and/or extend GR to overcome issues e.g., difficulties in quantizing GR and resolution of the curvature singularity, make it crucial to test the strong-field predictions of GR with the latest techniques and technologies.

Black holes (BH hereafter) are surprisingly ubiquitous objects in our Universe, and due to strong gravity regions in their neighborhoods, form the perfect candidates for testing theories of gravity. Within GR, under typical astrophysical conditions, a BH is an extremely simple object and its effect on the spacetime is described by very few parameters. Most commonly, these are the BH mass and spin and the object is known as a Kerr BH Kerr (1963). The assumption that astrophysical BHs are described by the Kerr metric is known as the Kerr hypothesis. (For the specific conditions and assumptions, see Ref. Chrusciel et al. (2012).) Alternative theories of gravity often introduce additional parameters, deforming the BH away from the Kerr solution. Observations of effects of BH have been a celebrated exercise in physics, given the potential for discovery of interesting phenomena. Some of the ways these observations are done are: X-ray spectroscopy (first measurements of the BH spin Dabrowski et al. (1997); Young et al. (1998)), gravitational wave interferometry (first observation of coalescence of a pair of BHs Abbott et al. (2016)), pulsar timing (first indirect detection of gravitational waves Taylor and Weisberg (1989)) and BH imaging (capturing an image of the region close to the BH horizon for the first time Akiyama et al. (2019b)).

In this work, our focus is on the technique of X-ray spectroscopy. In particular, we are interested in the reflection spectrum of BHs with accretion disks, which is in the X-ray band. Since the gravity of the BH affects its neighborhood and photons that travel from the neighborhood to us, the analysis of the reflection spectrum can be used to study the nature of BH itself Bambi (2017a); Bambi et al. (2018). The most advanced model for calculation of the reflection spectrum in the Kerr case is relxill García et al. (2014); Dauser et al. (2014). We have extended this model to relxill_nk, which can calculate the reflection spectrum for non-Kerr metrics Bambi et al. (2017); Abdikamalov et al. (2019a). The model has been applied to X-ray observations of several astrophysical BHs to place constraints on deviations away from the Kerr solution Cao et al. (2018); Tripathi et al. (2018); Wang et al. (2020); Xu et al. (2018); Choudhury et al. (2019); Zhou et al. (2018, 2019); Tripathi et al. (2019a, b); Zhang et al. (2019a); Tripathi et al. (2019c).

The non-Kerr metric most commonly used with relxill_nk is the one proposed in Ref. Johannsen (2013) (See also Vigeland et al. (2011)). This metric, referred to as the Johannsen metric hereafter, albeit not a solution of a well-defined alternative theory of gravity, preserves the symmetry of the Kerr metric associated with a Killing tensor Carter (1968); Papadopoulos and Kokkotas (2018). It can be mapped to BH solutions from some alternative theories with suitable combinations of the new parameters introduced in the metric. In this sense, it is a good candidate for theory-agnostic tests of the Kerr hypothesis. But it suffers from some important limitations. BHs in alternative theories do not generally have the additional symmetry of the Kerr solution. Moreover, the parametrization used in Ref. Johannsen (2013) for the deformation functions in the metric is based on an $M/r$ expansion (in units of $c=G=1$), where $M$ is the BH mass and $r$ the radial coordinate. This expansion becomes problematic at small $r$ for rapidly rotating BHs (e.g., near the inner edge of an accretion disk). This is exactly the regime where X-ray spectroscopy is most powerful Dauser et al. (2013); Abdikamalov et al. (2019b) and typically used Reynolds (2014); Nampalliwar and Bambi (2018). It is, therefore, essential to use a better metric to test the Kerr hypothesis, especially with X-ray reflection spectroscopy. In this paper, we present the application of a different metric to relxill_nk, which overcomes these limitations. This metric has been proposed in Ref. Konoplya et al. (2016), and we will refer to it as the KRZ metric hereafter.

The paper is organized as follows: in Sec. II, we review the KRZ metric and specify the deformation parameters, and review the relxill_nk framework and the numerical methods involved. Sec. III presents a qualitative analysis of the effect of the deformation parameters on relevant quantities (e.g., the innermost stable circular orbit) and observables (e.g., the reflection spectrum). In Sec. IV, the new model is applied to X-ray observations of a supermassive BH, Ark 564. The findings of the Ark 564 data analysis are discussed and compared with estimates from other observational techniques in Sec. V. Also discussed in Sec. V are issues with systematic uncertainties and future works. Throughout the paper, we will express quantities in units of $c=G=1$. Except when mentioned explicitly, we will set the BH mass $M=1$.

II Review

Non-Kerr metrics can be classified in two categories. The top-down metrics are those which are obtained as a solution of an alternative theory of gravity, e.g., the Einstein-dilaton-Gauss-Bonnet BHs Kanti et al. (1996); Kleihaus et al. (2016); Ayzenberg and Yunes (2014); Maselli et al. (2015); Kokkotas et al. (2017); Nampalliwar et al. (2018), the Chern-Simons BHs Yunes and Pretorius (2009); Yagi et al. (2012); McNees et al. (2016); Delsate et al. (2018), and the Kerr-Sen BHs Sen (1992); Hioki and Miyamoto (2008); Dastan et al. (2016); Uniyal et al. (2018). The bottom-up metrics on the other hand are obtained not from a specific theory of gravity but by generalizing the Kerr metric Vigeland et al. (2011); Johannsen (2013); Konoplya et al. (2016); Papadopoulos and Kokkotas (2018); Carson and Yagi (2020). Each approach has its advantages and disadvantages. Top-down metrics are difficult to obtain and might be known only numerically, but testing for them amounts to testing an alternative theory of gravity. Bottom-up metrics may have pathologies in the spacetime, but they can be mapped to several top-down metrics and thus constraints on parameters of bottom-up metrics translates to constraints on several top-down metrics. Our focus in the present work is on a particularly attractive bottom-up case, the KRZ metric.

II.1 The KRZ metric

The KRZ metric Konoplya et al. (2016) is based on a generic stationary and axisymmetric metric written in the Boyer-Lindquist like $(t,r,\theta,\phi)$ coordinates, where $t$ and $\phi$ are along the direction of a timelike and a spacelike Killing vector, respectively, $r$ and $\theta$ are orthogonal to each other and to $t$ and $\phi$, and at spatial infinity $(r,\theta,\phi)$ reduce to the standard spherical coordinates. The metric functions are written in terms of continued fraction expansions in the polar and radial coordinates.

The KRZ metric has several advantages over other bottom-up metrics (See Konoplya et al. (2016); Ni et al. (2016) for more details):

Parametrizations based on expansion in $M/r$ suffer from the inherent weakness that close to the horizon, where effect of gravity are strongest, higher-order parameters become equally important, making it impossible to isolate the dominant terms in the expansion. The continued fraction expansion trick allows for quicker convergence with fewer parameters. 2. 2.

A crucial feature of any bottom-up metric is its ability to map to several top-down metrics. The KRZ metric is generic enough to allow mapping to various top-down metrics with fewer parameters compared to other bottom-up metrics Kokkotas et al. (2017); Nampalliwar et al. (2018); Konoplya and Zhidenko (2020). 3. 3.

Several bottom-up metrics are based on the Janis-Newman transformation, which is not guaranteed to work beyond general relativity Hansen and Yunes (2013). The KRZ metric uses a different approach to avoid this.

The line element of the KRZ metric reads Konoplya et al. (2016); Ni et al. (2016)

[TABLE]

where each function $N,B,W$ and $K$ depends on $r$ and $\theta$, and is written in a way that separates the terms that are constrained asymptotically (e.g., through parametrized post-Newtonian constraints) and those that are constrained near the BH horizon. For instance,

[TABLE]

where $x=1-r_{0}/r$ and $r_{0}$ is the BH horizon radius in the equatorial plane, $y=\cos\theta$, and the $b_{i0}$ terms are constrained to ensure $B$ has the correct asymptotic behavior. The $\widetilde{B}_{i}(x)$ terms encode all the strong-field deviations of the metric and are written as Konoplya et al. (2016)

[TABLE]

Since we are interested in parametric deformations away from the Kerr solution in the strong field regime, we will set all asymptotically constrained terms to their GR values and only retain the terms calculated close to the BH horizon. We will follow the choice made in Ref. Ni et al. (2016) for the strong-field parameters to be retained. Assuming reflection symmetry across the equatorial plane and neglecting coefficients of higher orders, we get Ni et al. (2016)

[TABLE]

Here $a_{*}=J/M^{2}$ is the dimensionless spin parameter,

[TABLE]

and $\{\delta_{i}\}$ ($i=1,2,...6$) are six deformation parameters, which have been introduced to study deviations away from the Kerr metric. Note that the above expressions for the metric components are slightly different than those in Ref. Konoplya et al. (2016), since the expression given in Ref. Konoplya et al. (2016) do not reduce to the Kerr metric in the limit of zero deformation. The expressions above do reproduce the Kerr solution exactly when all $\delta_{i}$ are set to zero. The physical interpretation of the deformation parameters can be summarized as follows (see Ref. Konoplya et al. (2016) for more details):

[TABLE]

With this choice, the mass-quadrupole moment is the same as in the Kerr metric, and deviations from the Kerr solution are only possible in the strong gravity region.

Arbitrary deviations away from Kerr can create pathologies (e.g., closed time-like curves, Lorentz signature violation) in the spacetime. To avoid such scenarios, we impose restrictions on the above deformation parameters so that the following conditions are satisfied everywhere outside the horizon:

The metric determinant is always negative. 2. 2.

The metric coefficient $g_{\phi\phi}$ is greater than zero, and 3. 3.

$N^{2}$ is non-vanishing.

The explicit restrictions on each deformation parameter can be derived from these conditions using Eq. 1 and 4. Assuming only one non-zero $\delta_{i}$ at a time, we get the following restrictions on each $\delta_{i}$:

[TABLE]

Note that currently relxill_nk allows for variation in one deformation parameter at a time, so we will analyze the effect of each $\delta_{i}$ while fixing all other $\delta$s to zero. While mapping to BHs of alternative theories of gravity is unlikely to result in a single non-zero deformation parameter (e.g., see the mappings to static Einstein-dilaton-Gauss-Bonnet BHs in Ref. Kokkotas et al. (2017) and slowly rotating Einstein-dilaton-Gauss-Bonnet BHs in Ref. Konoplya et al. (2016)), we consider our analysis as a test of the validity of GR, in the spirit of null tests Psaltis et al. (2015); Ashby et al. (2018).

II.2 The Xspec model

The standard astrophysical system in X-ray reflection spectroscopy is modeled with the BH-disk-corona model Fabian et al. (1989); Laor (1991); Bambi (2017b); Nampalliwar and Bambi (2018). The BH is assumed to be surrounded by a geometrically thin and optically thick disk Novikov and Thorne (1973),111Various other disk structures are possible. Studies of their effect on the reflection spectrum in the presence of non-Kerr metrics is an ongoing effort Riaz et al. (2020a); Abdikamalov et al. (2020). with its inner edge at some radius $r_{\textrm{in}}$, often at the innermost stable circular orbit (ISCO), and the outer edge at some large radius $r_{\textrm{out}}$. In addition, the system possesses a “corona”, which is thought to be a cloud of hotter (relative to the disk) gas or the base of an astrophysical jet. Fig. 1 illustrates the system under discussion. The disk emits as a blackbody locally and as a multi-temperature blackbody when integrated radially (labeled the thermal component in Fig. 1). Inverse Compton scattering of the thermal component by the corona produces X-rays (labeled the power-law component in Fig. 1), some of which returns to the disk and is reflected (labeled the reflected component in Fig. 1). Among these, the reflected component is most promising for studying the effects of strong gravity Reynolds (2014); Bambi (2017a).

Modeling the reflection component requires some understanding of the various physical parameters of disk-corona model. relxill_nk has several parameters to account for the different aspects of the system. These include the inner and outer edge of the disk, inclination of the disk relative to the observer, the disk’s elemental constitution and their ionization, and the emissivity profile. The emissivity profile determines the reflection spectrum at the source, i.e., at the disk, and depends strongly on the coronal geometry. Since the latter is poorly understood, unless in specific cases like the lamp-post geometry, the emissivity profile is modeled by a power-law (intensity $\propto 1/r^{q}$) or a broken power-law (intensity $\propto 1/r^{a}$ for $r\leq r_{\textrm{br}}$ and $\propto 1/r^{b}$ for $r>r_{\textrm{br}}$). Apart from the BH spin, the spacetime is determined by two parameters associated with deviations from Kerr: one parameter picks a specific deformation parameter, e.g., KRZ $\delta_{1}$; the other parameter defines the value of the chosen deformation parameter.222More details can be found on the public version webpage at Refs. rel (a, b). This flexibility allows us to the use the same framework of relxill_nk to analyze deformation parameters of various metrics. Note that relxill_nk allows one deformation parameter at a time to vary, a limitation imposed by the FITS table file (see below) which grows exponentially with each additional deformation parameter.

II.3 Numerical method

The numerical method has been presented in detail earlier Bambi et al. (2017); Abdikamalov et al. (2019a). Here we briefly review it. The relxill_nk framework uses the transfer function formalism, introduced in  Cunningham (1975) and implemented in the relxill suite of models Dauser et al. (2013); Garcia et al. (2013); García et al. (2014) to optimize the calculation of the X-ray flux. In terms of the transfer function, the flux can be written as Bambi et al. (2017):

[TABLE]

Here $I_{o}$ and $I_{e}$ are the specific intensities of the radiation detected by the distant observer and the emitter respectively. These are related via the Liouville’s theorem: $I_{o}=g^{3}I_{e}$, where $g=\nu_{o}/\nu_{e}$ is the redshift factor, $\nu_{o}$ is the photon frequency as measured by the distant observer, and $\nu_{e}$ is the photon frequency in the rest frame of the emitter. $r_{e}$ is the emission radius in the disk and $\vartheta_{e}$ is the photon’s direction relative to the disk at the point of emission. $X$ and $Y$ are the Cartesian coordinates of the image of the disk in the plane of the distant observer, $D$ is the distance of the observer from the source, and $d\tilde{\Omega}=dXdY/D^{2}$ is the element of the solid angle subtended by the image of the disk in the observer’s sky. The transfer function itself is defined as

[TABLE]

Here $g^{*}$, the normalized redshift factor, is defined as

[TABLE]

where $g_{\rm max}=g_{\rm max}(r_{e},i)$ and $g_{\rm min}=g_{\rm min}(r_{e},i)$ are, respectively, the maximum and the minimum values of the redshift factor $g$ at a constant $r_{e}$ and for a given viewing angle of the observer. The $r_{e}$-integral ranges from the inner to the outer edge of the disk, whereas the $g^{*}$-integral ranges from 0 to 1. Since there are two ways to go from $g^{*}=0$ to $g^{*}=1$ along a constant $r_{e}$ ring, there are two branches of the transfer function and the integral in Eq. 17 needs to be performed along both the branches. Furthermore, if the emission from the disk is not isotropic, then we also need to calculate $\vartheta_{e}$ at each $r_{e}$ and $g^{*}$ along each branch, and Eq. 17 will be

[TABLE]

The transfer function separates the spacetime effects like the motion of gas and photons (encoded in the transfer function) from the local microphysics (encoded in the specific intensity at the emission point), and thus acts as an integration kernel for the calculation of flux. Such a separation enables quick computation of the reflection spectrum from a grid of transfer functions for any intensity profile, without the need to retrace photon trajectories, making analysis of X-ray reflection data with relxill_nk possible in practice.

These transfer functions and emission angles are computed once and stored in a FITS (Flexible Image Transport System) table, and only the integral in Eq. II.3 needs to be calculated during data analysis. The table has a grid in three dimensions: spin $a_{*}$ (30 values), deformation parameter $\delta_{i}$ (30 values), and inclination angle $i$ (22 values). The grid points in the $a_{*}$ and $i$ dimension are non-uniform, mutually independent and follow the scheme in relxill. In particular, $a_{*}$ ranges from $-0.998$ to $0.9982$ and $i$ ranges from $3^{\circ}$ to $89^{\circ}$. The range in the $\delta_{i}$ dimension, in some cases, depends on the spin parameter: this is to ensure that the restrictions on $\delta_{i}$, given in Sec. II.1, are followed. We use the following scheme for all $\delta_{i}$: The largest extent of $\delta_{i}$ at each spin is $-5$ to $5$, distributed uniformly between the two bounds. If the restriction given in Sec. II.1 is stronger than this bound, the range of that $\delta_{i}$ at that spin is set according to the restriction. Examples of this scheme can be seen in Fig. 2, where the disallowed region is marked with stripes.

At each grid point (namely a specific ($a_{*}$, $\delta_{i}$, $i$) on the grid described above), the accretion disk is discretized in 100 values of $r_{e}$, between $r_{\rm ISCO}$ and $1000M$, and 20 values of $g^{*}$ on each branch, between $10^{-3}$ and $1-10^{-3}$.333Since the Jacobian in the transfer function diverges at $g^{*}=0$ and $1$, we offset the values by $10^{-3}$. Photons are back-traced from the observer plane (placed at a large distance where the spacetime is effectively flat) to the accretion disk by solving the geodesic equations using a fourth-order Runge-Kutta scheme. An adaptive algorithm fine-tunes the coordinates on the observer plane so that the photon when back-traced lands at the exact $r_{\textrm{e}}$. For each such “central” photon, the redshift and emission angle are calculated, and four photons, closely spaced in the observer plane, are launched to calculate the Jacobian and subsequently the transfer function. An interpolation routine then computes these quantities at the 20 values of $g^{*}$ and stores them in the FITS table.

III Qualitative analysis of deviation parameters

We now look at the qualitative effect of the deformation parameters on some typical observables and quantities in X-ray reflection spectroscopy, to get some intuition how the observables are affected by the different kind of deformations, and anticipate when, if at all, they can be used to estimate the parameters.

One of the most important quantities is the ISCO. Since the emissivity profile has a radial fall-off, radiation from the innermost region has a significant influence on the X-ray spectra. The influence of the deformation parameters drops off radially as well, manifesting most prominently close to the horizon. Theoretically, the inner edge of the disk is assumed to be at ISCO within the Novikov-Thorne thin-disk models Novikov and Thorne (1973). Observationally too, disks in thermal state and accreting at $\sim 5\%-30\%$ of the Eddington limit are expected to be bounded by the ISCO McClintock et al. (2006); Steiner et al. (2010). Thus, the ISCO serves as a proxy for determining constraints on deformation parameters with typical X-ray data. For axisymmetric metrics, the calculation of the ISCO radius involves only the $g_{tt}$, $g_{\phi\phi}$ and $g_{t\phi}$ components of the metric in the equatorial plane (see, for instance, Appendix A.3 in Bambi et al. (2017)). Thus, among the six deformation parameters under consideration here, only $\delta_{1}$ and $\delta_{2}$ influence the ISCO radius (see Eqs. 1), while $\delta_{3}$, $\delta_{4}$, $\delta_{5}$, and $\delta_{6}$ have no effect (but they do affect other quantities like the redshift and the reflection spectrum, as we will see shortly). Fig. 2 shows how the ISCO radius changes with $\delta_{1}$ and $\delta_{2}$ in the top and bottom panels, respectively. Since the BH spin also has a strong effect on the ISCO, and in X-ray spectroscopy we typically measure both spin and deformation, each panel shows the variation in the ISCO radius on a spin-deformation parameter axes ($a_{*}-\delta_{1}$ in the top panel, $a_{*}-\delta_{2}$ in the bottom panel). In both cases, small ISCOs ($r_{\textrm{ISCO}}\lesssim 2M$) are confined to high spins, implying good estimates of $a_{*}$ are possible with X-ray data. For larger ISCO, the contours are spread over large ranges of $a_{*}$, making good spin estimates unlikely. For both $\delta_{1}$ and $\delta_{2}$, the ISCO contours stretch over a wide range, making precise estimation very difficult. The only exception is very small ISCOs ($r_{\textrm{ISCO}}\lesssim 1.5M$) for which both $\delta_{1}$ and $\delta_{2}$ are confined to a small range.

We now look at the primary observable in X-ray reflection spectrocopy and the output of the relxill_nk model. Fig. 3 shows the reflection component of the flux from a typical BH-disk system, for various $\delta_{i}$s. Each panel shows one specific $\delta_{i}$, and within each panel spectra for three different values of the $\delta_{i}$ parameter are shown. Remaining model parameters are kept constant in all panels (and listed in the figure caption). In particular, $a_{*}$ is set to $0.99$ and $q$ to $6$, since X-ray reflection spectroscopy is most useful for performing tests of gravity in systems with very high BH spins and steep emissivity profiles Fabian et al. (2002); Dauser et al. (2013). In each panel, we plot the spectra for zero deformation ($\delta_{i}=0$) and positive/negative deformation ($\delta_{i}\lessgtr 0$). The positive/negative deformation value is typically set to $\pm 2.5$, except when the range of $\delta_{i}$ restricts it to a smaller deformation (see, e.g., Fig. 2 and Sec. II.3). Several features can be seen in the plots. The spectrum is very sensitive to positive $\delta_{1}$ while negative values of $\delta_{1}$ are severely restricted and the spectrum does not change much for the most negative $\delta_{1}$ allowed. In the $\delta_{2}$ case, the spectrum is very sensitive to negative values of $\delta_{2}$, while positive values have a moderate effect. $\delta_{3}$ also has a moderate-to-strong effect on the spectrum for negative values, whereas positive $\delta_{3}$ affects the spectrum rather weakly. $\delta_{4}$ has negligible effect for either positive or negative values, but $\delta_{5}$ has a strong effect for positive values. Finally, $\delta_{6}$ has a strong effect on the spectrum for negative values, whereas positive values do not affect the spectrum much. In each panel, we also plot the spectra from relline Dauser et al. (2010), which includes only the Kerr metric, to facilitate comparison of the zero deformation case with the Kerr case. For all $\delta_{i}$s, the zero deformation case matches very well with the Kerr case.

The ability of X-ray spectroscopy, and in fact any technique, to constrain parameters depends not only on the effect of the parameter on the relevant observable, but also on the degeneracy among the parameters in terms of their effect on the observable. In particular, above we looked at the effect of $\delta_{i}$ on the reflection spectra at fixed BH spin. But it is well known that spin and deformation parameter affect the spectra in a somewhat similar way, resulting in a degeneracy. For statistically robust parameter estimation, it is critical to ensure all degeneracies are accounted for. We will see instances of this in the next sections where we estimate $\delta_{i}$ from astrophysical data.

IV Spectral analysis of Ark 564 data

In this section, we use our newly developed model to analyze an X-ray observation of the supermassive BH in Ark 564 Walton et al. (2013).

IV.1 Review

Ark 564 is classified as a narrow line Seyfert 1 galaxy at redshift $z=0.0247$. Since first observations with XMM-Newton in 2000/2001 Vignali et al. (2004), it has been studied by several authors since it appears as a very bright source in the X-ray band Papadakis et al. (2007); McHardy et al. (2007); Kara et al. (2017); Barua et al. (2020). It is a good first candidate for tests of general relativity with a new metric for the following reasons. Firstly, previous studies have shown that the inner edge of the disk may be very close to the central object, which maximizes the signatures of the strong gravity region Walton et al. (2013). Secondly, the source has a simple spectrum. There is no obvious intrinsic absorption to complicate the determination of the reflected emission. The same dataset was analyzed by some of us in Ref. Tripathi et al. (2018) to constrain the Johannsen parameters $\alpha_{13}$ and $\alpha_{22}$, and in Ref. Tripathi et al. (2019c) to constrain the Johannsen parameter $\epsilon_{3}$. Therein, this source and this dataset was found to be easy to analyze and provided good constraints on the deformation parameters. There are other X-ray sources with the potential to provide stronger constraints than this data (e.g. GRS 1915+105 Zhang et al. (2019b), MCG-06-30-15 Tripathi et al. (2019a)) and their analysis with the KRZ metric is underway.

IV.2 Observations and data reduction

Suzaku observed Ark 564 on 26-28 June 2007 (Obs. ID 702117010) for about $80$ ks. For low energies ($<10$ keV), Suzaku has four co-aligned telescopes which are used to collect photons onto its CCD detectors X-ray Imaging Spectrometer (XIS). XIS is comprised of four detectors; XIS0, XIS2, and XIS3 are front-illuminated and XIS1 is back-illuminated. We only used data from the front-illuminated chips because XIS1 has a lower effective area at 6 keV and a higher background at higher energies. XIS2 data were not used in our analysis because of the anomaly after 9 November 2006.

We have used HEASOFT version 6.24 and CALDB version 20180312 for the data reduction. AEPIPELINE script of the HEASOFT package has been used for reprocessing and screening of the raw data. The ftool XSELECT and ds9 were used to extract the XIS source and background spectrum from a 3.5 arc minute radius. The background region was selected as far as possible from the source so as to avoid contamination. The RMF file was generated using XISRMFGEN and the ARF file which corresponds to effective area of the telescope was generated using XISSIMARFGEN. At last, the final source spectrum, background spectrum and response file were generated by combining the data from XIS0 and XIS3 using ADDASCASPEC. The data was then grouped using GRPPHA to get minimum 50 counts per bin so as to use $\chi^{2}$ statistics in our spectral analysis. We have also excluded the energy range between 1.7-2.5 because of calibration issues.

IV.3 Modelling and results

We fit the data with Xspec v12.10.0c. The data has been analyzed before, with models that incorporate the Kerr metric Kara et al. (2017) and with those incorporate non-Kerr metrics Tripathi et al. (2018, 2019c). In particular, we will use the findings of Ref. Tripathi et al. (2018), where a comprehensive analysis has been done of the model combinations, and the optimal combination found had a relativistic and a non-relativistic reflection component. Of course, there is no a priori reason to assume that the models used to fit some data must be the same for all non-Kerr parameters, so we follow the standard approach by first fitting the power-law component (see Sec. II.2 and Fig. 1), followed by the reflection components.

Model 1

tbabs*zpowerlw

tbabs describes the galactic absorption Wilms et al. (2000) and we fix the galactic column density to $N_{H}=6.74\cdot 10^{20}\textrm{cm}^{-2}$ gal . zpowerlw describes a redshifted photon powerlaw spectrum. The best-fit values are reported in the second column in Tab. 1, and the data-to-model ratio plots are presented in the top panel in Fig. 4. An excess of photons at low energies and a broad iron line around 6.4 keV can be seen in the panel, suggesting the presence of a reflected component. Invoking this paradigm, we proceed to fit the spectrum using the reflection models relxill_nk and xillver.

Model 2

tbabs*(zpowerlw + relxill_nk + xillver)

The power-law component is still modeled with zpowerlw as before, while the reflection is modeled with relxill_nk and xillver. For the fit, the disk is assumed to extend from the ISCO to $400M$. The emissivity profile is assumed to be a simple power law, whose index $q$ is a free parameter. The photon index in relxill_nk is tied to zpowerlw and the cut-off energy frozen to the default value ($300$ keV). The ionization and iron abundance of the disk, along with the BH spin and deformation parameter, are kept free. The parameters of xillver are mostly fixed or tied (ionization fixed to zero, iron abundance, cut-off energy, inclination tied to relxill_nk, photon index tied to zpowerlw), the only free parameter being the norm. The fit is performed six times, once for each $\delta_{i}$.

During the fit, inclination had to be handled with care, for the following reason: previous constraints on inclination with this (or combined with this) observation of Ark 564 have not been very good, either unbounded from below (see Tab. I of Ref. Tripathi et al. (2018) and Tab. 1 of Ref. Kara et al. (2017)) or very precise but inconsistent (see Tab. II of Ref. Tripathi et al. (2018), Tab. II of Ref. Tripathi et al. (2019c) and Tab. 1 of Ref. Kara et al. (2017)). Keeping inclination free in our analysis, we run into similar issues, i.e., either the parameter is largely unbounded or stuck at the smallest allowed value ($\sim 3^{\circ}$) with little uncertainty. Instead, we bound it to lie around $30^{\circ}$ with a $10\%$ uncertainty (i.e., $27^{\circ}<i<33^{\circ}$), which encompasses to a large extent estimates obtained with a different observation of Ark 564 Kara et al. (2017); Barua et al. (2020) made with NuSTAR, and with the Suzaku observation and different deformation parameters Tripathi et al. (2018, 2019c).

The data-to-model ratios are plotted for each iteration of Model 2, i.e., for each $\delta_{i}$, in Fig. 4. For each best-fit model, parameter values are listed in Tab. 1, with associated uncertainty at $90\%$ confidence level (CL). We find that, in all cases, the fit is good with a $\chi_{\textrm{red}}^{2}$ close to 1 and no major unmodeled features in the ratio plots. We find a high value for the photon index $q$, that is, most of the radiation seems to come from very inner part of the accretion disk, consistent with previous analyses Tripathi et al. (2018, 2019c). The photon index is slightly higher than Ref. Kara et al. (2017) (both being higher than in Ref. Barua et al. (2020)) but is consistent for each $\delta_{i}$ case. Ionization is consistent amongst the Model 2 iterations and close to previous analyses (lower than Ref. Tripathi et al. (2019c), higher than Ref. Tripathi et al. (2018)), The spin parameter $a_{*}$ is always very close to 1. This is consistent with the previous analysis results of Ark 564 and suggests that it is a high spin BH. Iron abundance is close to solar and matches well with Ref. Tripathi et al. (2018) (see also Ref. Kara et al. (2017)). The inclination is, in most cases, pegged to the smallest allowed value ($27^{\circ}$) and, in some cases, unbounded within its restricted range.

V Discussion

The primary aim of this analysis is to study the spacetime metric around the supermassive BH in Ark 564 and estimate the deformation parameters $\delta_{i}$. Tab. 1 lists the best-fit value of each $\delta_{i}$ and the associated uncertainty (at $90\%$ CL). In particular,

[TABLE]

As we saw in the qualitative discussion of Sec. III, different parameters affect the reflection spectrum differently and thus, using the reflection data of Ark 564, are constrained to different levels. Specifically, $\delta_{1}$ is strongly constrained on the positive side whereas on the negative side the bound comes from the spacetime regularity conditions of Sec. II.1. $\delta_{2}$ is strongly constrained on the negative side (considering the best-fit value of $\delta_{2}$, $-0.25$) and moderately constrained on the positive side, while $\delta_{3}$ is constrained strongly on the negative side and almost unconstrained on the positive side. $\delta_{4}$ has almost no effect on the reflection spectrum and consequently is very poorly constrained, in fact the upper bound exceeds the grid size which is at $\delta_{4}=5$. The best-fits for $\delta_{5}$ and $\delta_{6}$ are both quite far from [math], though both have large uncertainties with $\delta_{5}$ hitting the theoretical limit (based on regularity of spacetime) and $\delta_{6}$ very close to the computational limit (at $\delta_{6}=-5$). In both $\delta_{5}$ and $\delta_{6}$ cases, the Kerr solution is recovered at $99\%$ CL (shown in parentheses).

The KRZ metric has been used to perform tests of GR with other experimental techniques as well. Gravitational waves, in particular, have been used extensively in recent years to test GR in the strong field regime. The spherically symmetric version of the KRZ metric Rezzolla and Zhidenko (2014) has been used in Ref. Cardenas-Avendano et al. (2020a) to put constraints on one of its deformation parameters, using data from the first two observing runs of LIGO/VIRGO. The parameter analyzed in Ref. Cardenas-Avendano et al. (2020a), $a_{1}$, is the same as $\delta_{1}$ of the present work. Since all the other parameters in that work are set to their GR value, the estimates obtained there can be faithfully compared with those obtained in the present work. Comparing the constraint on $\delta_{1}$ from above with Tab. 1 in Ref. Cardenas-Avendano et al. (2020a), we see that the constraints of the present work are of the same order as Ref. Cardenas-Avendano et al. (2020a): better than those obtained from the gravitational wave event GW170104 and worse than the other four events. It is important to remember though that the study in Ref. Cardenas-Avendano et al. (2020a) was done under the assumption of spherically symmetric BHs, which certainly will not be valid for all the gravitational wave events, and inclusion of spin will worsen the estimates there. Another recent work Psaltis et al. (2020) uses the results from the observation of the supermassive BH at the center of the M$87$ galaxy, reported by the Event Horizon Telescope collaboration Akiyama et al. (2019b), to place constraints on the leading order strong-gravity deformation parameter, same as our $\delta_{1}$. The constraints obtained there are a few times weaker than those reported here. The constraints are summarized in Tab. 2.

It is well known that model parameters can be degenerate, so while estimating them one needs to be careful in handling these degeneracies. We take care of this by marginalizing over all the free model parameters, and the uncertainties reported in Tab. 1 are therefore robust statistically (systematic uncertainties are a different story, see the discussion later). In particular, there exists a strong degeneracy between $a_{*}$ and the deformation parameters, since both appear in the metric and, consequently, affect the spacetime in a somewhat similar way. Therefore, any test of GR has to contend with the fact that simultaneous measurement of both spin and non-Kerr parameters will have an intrinsic degeneracy, resulting in weaker constraints on each parameter (than if the BH was not spinning or if the non-Kerr parameters were identically zero). Fig. 5 shows this for the fits made with the Ark 564 data discussed above444Specifically, after obtaining the best-fit model in each $\delta_{i}$ case, we use Xspec’s steppar command over two dimensions, $a_{*}$ and $\delta_{i}$, which steps through a series of pairs of ($a_{*},\delta_{i}$) values, finds a new fit at each step by marginalizing over all other free parameters, and reports the $\chi^{2}$ at each step. The contours in the plots represent the change in $\chi^{2}$ at each step relative to the $\chi^{2}$ of the original best-fit., with each panel showing the relative degeneracy between $a_{*}$ and one of the $\delta_{i}$. The red, green, and blue lines indicate, respectively, the 68%, 90%, and 99% confidence level contours. Generically, we can say that the constraints on each $\delta_{i}$ would be very different if the degeneracy with spin was not taken into account. For instance, we see in the panel with $\delta_{6}$ (bottom right panel) that at fixed spin ($a_{*}=0.99$), $\delta_{6}=-2$ is distinct from $\delta_{6}=0$, excluded at more than $99\%$ CL. This agrees with the qualitative analysis of Sec. III. But the analysis there was limited to constant BH spin. When the spin is allowed to vary, a larger range of $\delta_{6}$ ($\delta_{6}=-4$ at $a_{*}=0.97$) becomes degenerate. The degeneracy manifests in another way in the contours between $a_{*}$ and $\delta_{1}$ (top left panel) and $a_{*}$ and $\delta_{2}$ (top right panel), respectively. As discussed in Sec. III, the ISCO contour can be used as a proxy for the degeneracy, especially near the Kerr solution. In both the panels, we find that the degeneracy contours follow the constant ISCO contours close to the Kerr solution (see Fig. 2). For larger deviations away from Kerr, the spectrum (which is affected by several things and not just the ISCO location) becomes distinct, and the ISCO ceases to be the driver of degeneracy.

One of the aspects beyond the scope of this work that is important for parameter estimation is systematic uncertainties. They arise from the fact that the analysis makes a series of assumptions. The model, for example, assumes a thin disk described by the Novikov-Thorne model Novikov and Thorne (1973). For astrophysical systems, this is a valid assumption only when the accretion rate is not too high (nor too low) Penna et al. (2010); Kulkarni et al. (2011). Investigations with thick disks have shown that this assumption can potentially cause systematic error in spin up to $0.2$ Taylor and Reynolds (2018) and $0.4$ Riaz et al. (2020b), and more than $1$ for deformation parameters Riaz et al. (2020a). Another assumption, that particles move in circular orbits on the disk can also lead to systematic errors Tripathi et al. (2020). Higher order effects in the context of disk reflection, neglected in the current models, can be important in some cases Zhu et al. (2012); Zhou et al. (2020); Cardenas-Avendano et al. (2020b). The model is limited to only one non-zero deformation parameter. Degeneracy among the various deformation parameters is, therefore, beyond the scope of the model. While this is fine for verifying GR, mapping to BH metrics from alternative theories of gravity will certainly involve multiple deformation parameters. In this context, a recent study based on gravitational wave ringdown Völkel and Barausse (2020) has shown the severe weakening of constraints on each when multiple deformation parameters are set free.

We emphasize that this work is an initial exploration with several aspects to be explored in future. These include, but are not limited to, studies involving more parameters (e.g., estimating $\epsilon_{0}$ parameter, see Eq. 1), estimating multiple deformation parameters simultaneously (e.g., as done in Völkel and Barausse (2020)), employing more sophisticated disk geometry (e.g., as done in Abdikamalov et al. (2020)), studying more astrophysical sources (e.g., MCG-06-30-15 Fabian et al. (2002); Marinucci et al. (2014); Tripathi et al. (2019a), GRS 1915+105 McClintock et al. (2006); Miller et al. (2013); Zhang et al. (2019b), GX 339-4 Reis et al. (2008); García et al. (2015); Wang et al. (2020)), and exploring potential synergy with other observational techniques Cardenas-Avendano et al. (2020a); Nampalliwar et al. (2020); Völkel and Barausse (2020).

Acknowledgments – We thank Alejandro Cardenas-Avendano for pointing out the issue with matching the metric to the Kerr solution, and Roman Konoplya and Alexander Zhidenko for providing the correct metric components. We also thank the anonymous referee whose valuable feedback helped us improve the paper. S.N. acknowledges support from the Excellence Initiative at Eberhard-Karls Universität Tübingen. S.N. and J.A.G. also acknowledge support from the Alexander von Humboldt Foundation. A.B.A. acknowledges the support from the Shanghai Government Scholarship (SGS). A.T. acknowledges support from the China Scholarship Council (CSC), Grant No. 2016GXZR89. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG. This work was supported by the National Natural Science Foundation of China (NSFC), Grant No. U1531117, and Fudan University, Grant No. IDH1512060.