Testing the Kerr Black Hole Hypothesis Using X-Ray Reflection Spectroscopy and a Thin Disk Model with Finite Thickness

Einstein's theory of general relativity is a pillar of modern physics and in agreement with all the available observational tests (Will 2014). However, the theory has been primarily tested in weak gravitational fields, while its predictions in the strong field regime have only recently being put to test. Astrophysical black holes are ideal laboratories for testing general relativity in the strong field regime and a number of theoretical reasonings point to the possibility that the spacetime metric around these objects can present macroscopic deviations from the predictions of Einstein's gravity (see, for instance, Dvali & Gomez 2013; Giddings & Psaltis 2018; Carballo-Rubio et al. 2020).

In four-dimensional general relativity, uncharged black holes are relatively simple systems. They are described by the Kerr solution (Kerr 1963) and are completely specified by only two parameters, representing, respectively, the mass M and the spin angular momentum J of the black hole. This is the well-known conclusion of the no-hair theorems, and it holds under specific assumptions (Carter 1971; Robinson 1975; Chruściel et al. 2012). It is also quite remarkable that the spacetime metric around an astrophysical black hole formed from the complete collapse of a progenitor body should be well approximated by the simple Kerr solution. For example, the presence of an accretion disk or of a nearby star has a very small impact on the near horizon metric and can normally be ignored (Bambi et al. 2014; Bambi 2018). The search for possible deviations from the Kerr geometry in the strong gravity region of an astrophysical black hole can thus be a tool to constrain and find new physics.

The Kerr black hole hypothesis can be tested by studying the properties of the electromagnetic radiation emitted by material orbiting a black hole (Johannsen 2016; Bambi 2017; Krawczynski 2018; Zhou et al. 2019b). Among all the electromagnetic techniques for testing the near horizon region of black holes, X-ray reflection spectroscopy (Fabian et al. 1989; Brenneman & Reynolds 2006; Reynolds 2014) is the most mature one and the only one that can currently provide quantitative constraints on the black hole strong gravity region (see, for instance, Cao et al. 2018; Tripathi et al. 2019a, 2019b; Zhang et al. 2019a). Like any astrophysical measurement, even for X-ray reflection spectroscopy it is crucial to have a sufficiently sophisticated astrophysical model in order to limit the modeling systematic uncertainties.

X-ray reflection spectroscopy refers to the analysis of the features of the reflection spectrum of accretion disks. Our system is a central black hole accreting from a geometrically thin and optically thick disk, with the inner edge of the disk at the innermost stable circular orbit (ISCO). Similar disks are thought to form when the source is in the thermal state with an accretion luminosity between a few percent and about 30% of its Eddington limit (McClintock et al. 2006; Penna et al. 2010; Steiner et al. 2010). The gas of the accretion disk is in local thermal equilibrium and at any point on the surface of the disk the emission is like that of a blackbody. The spectrum of the whole disk is a multi-temperature blackbody-like spectrum because the temperature increases as the gas falls into the gravitational well of the black hole (Page & Thorne 1974; Zhang et al. 1997). The thermal emission of the accretion disk normally peaks in the soft X-ray band (0.1–1 keV) for stellar-mass black holes and in the optical/UV band (1–100 eV) for supermassive ones, as the disk temperature scales as M^−0.25 (Zhang et al. 1997). The "corona" is some hotter (∼100 keV), usually compact and optically thin, gas near the black hole. Thermal photons from the disk can inverse Compton scatter off free electrons in the corona, producing a power-law component with an exponential cutoff in the X-ray spectrum of the black hole (Sunyaev & Truemper 1979). The Comptonized photons can illuminate the accretion disk, producing the reflection component (George & Fabian 1991; Ross & Fabian 2005; García et al. 2013). The latter is characterized by fluorescent emission lines below 8 keV, notably the iron Kα complex at 6.4–6.79 keV, depending on the ionization of iron ions, and the so-called Compton hump peaking at 20–30 keV.

A relativistic reflection model relies on a model to calculate the reflection spectrum at every emission point on the disk (assuming Einstein's equivalence principle holds, these calculations only involve atomic physics) as well as on a disk-corona model⁶ and a spacetime metric, which are both necessary to calculate the reflection spectrum at the detection point far from the source. All these pieces have a number of parameters and a variation in the value of these model parameters can have an impact on the predicted reflection spectrum of an accreting black hole. Fitting observational data with the theoretical model, we can infer the value of the model parameters and thus the properties of the system. If we employ a spacetime metric with some parameters quantifying deviations from the Kerr spacetime, we can attempt to constrain possible deviations from the Kerr metric by fitting X-ray data of some reflection-dominated source with our model.

relxill_nk is a relativistic reflection model to test the Kerr black hole hypothesis (Bambi et al. 2017; Abdikamalov et al. 2019). It is an extension of the relxill package (Dauser et al. 2013; García et al. 2013, 2014) to non-Kerr spacetimes. As in relxill, in relxill_nk the reflection spectrum in the rest frame of the disk is modeled by xillver, the accreting matter is described by an infinitesimally thin Novikov–Thorne disk (Novikov & Thorne 1973; Page & Thorne 1974), and the disk's intensity profile is either described by a broken power law or is the profile generated by a corona with lamppost geometry. relxill_nk differs from relxill only in the spacetime metric. The main version of relxill_nk employs the Johannsen metric (Johannsen 2013), which is not an exact solution of some specific gravity model but a parametric black hole spacetime. The Johannsen metric has an infinite number of "deformation parameters" that quantify deviations from the Kerr background. With the spirit of a null experiment, we can fit the reflection spectrum of a source with relxill_nk, determine the values of the deformation parameters, and thus verify if they are consistent with the hypothesis that the metric around the source is described by the Kerr solution as required by general relativity. As it has been constructed, relxill_nk can easily employ any stationary, axisymmetric, and asymptotically flat metric in analytic form (see, e.g., Zhou et al. 2018; Nampalliwar et al. 2019; Tripathi et al. 2020; Zhou et al. 2020a).

Like in any astrophysical measurement, even for the tests of the Kerr metric with relxill_nk, it is crucial to limit the systematic uncertainties. Otherwise, in the presence of high-quality data, we could obtain precise but inaccurate measurements of the spacetime metric around an accreting black hole and our analysis may find deviations from the Kerr solution that, actually, are due to systematic uncertainties. Among all the systematic uncertainties, modeling uncertainties are normally the dominant ones. relxill_nk has a number of modeling uncertainties, ranging from simplifications in the nonrelativistic reflection model and in the disk-corona model to relativistic effects not taken into account (see, for instance, the discussion in Liu et al. 2019; Zhou et al. 2020b).

All the available relativistic reflection models assume that the black hole accretion disk is geometrically thin and that there is no emission of radiation inside the inner edge of the disk. For example, if we apply these models to sources accreting near their Eddington limit, the spin parameter can be easily overestimated (Riaz et al. 2020a, 2020b). Moreover, the accretion disk is always approximated as infinitesimally thin. For a real accretion disk, we should expect that the disk has a finite thickness and that the latter increases as the mass accretion rate increases. Employing a model with an infinitesimally thin accretion disk inevitably leads to modeling bias in the final measurements of some model parameters. The impact of such a simplification has been ignored for a long time and only recently Taylor & Reynolds (2018a) have presented a relativistic reflection model in which the accretion disk has a finite thickness. In the present paper, we implement the accretion disk model of Taylor & Reynolds (2018a) into rellxil_nk as a step of our program of developing relativistic reflection models for testing the Kerr black hole hypothesis in order to try to create a tool for precision tests of general relativity in the strong field regime. The implementation of a disk of finite thickness in rellxil_nk can be useful to analyze those sources with thicker accretion disks and providing more precise measurements of the deformation parameters of the spacetime.

The paper is organized as follows. In Section 2, we review the accretion disk geometry proposed in Taylor & Reynolds (2018a) and, in Section 3, we employ such a disk geometry in relxill_nk. In Section 4, we use the new model with a disk of finite thickness to analyze a Suzaku observation of the X-ray binary GRS 1915+105 and to explore the impact of the disk thickness in our tests of the Kerr metric. This Suzaku observation of GRS 1915+105 was studied in Zhang et al. (2019a) and currently provides one of the most precise measurements of the deformation parameters of the spacetime, so it is presumably quite sensitive to systematic uncertainties. Summary and conclusions are reported in Section 5. In the Appendix, we briefly review the Johannsen metric and its black hole parameter space. Throughout the paper, we use units in which G_N = c = 1 and a metric with signature (−+++).

We consider the accretion disk geometry proposed in Taylor & Reynolds (2018a). We assume that the accretion disk midplane lies in the θ = π/2 plane and is a radiation-pressure dominated, geometrically thin, and optically thick disk with a pressure scale height H defined as (Shakura & Sunyaev 1973)

$$\begin{eqnarray}&&H=\displaystyle \frac{3}{2}\displaystyle \frac{1}{\eta }\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}\right)\left[1-\sqrt{\displaystyle \frac{{r}_{\mathrm{ISCO}}}{\rho }}\right],\end{eqnarray} \tag{ 1 }$$

where ρ = r sin θ is the pseudo-cylindrical radius, $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ is the Eddington-scaled mass accretion rate, and r_ISCO is the ISCO radius, which is also the inner edge of the disk. The radiative efficiency is η = 1 − E_ISCO, where E_ISCO is the specific energy of a test particle in the midplane at r_ISCO. We assume that the surface of the disk is determined by the half thickness $z(\rho )=2H$ and that the disk rotates cylindrically ($\dot{\theta }=0$), which means that all matter at some pseudo-cylindrical radius ρ in the disk will have the same orbital velocity as the material at the same cylindrical radius in the equatorial plane (θ = π/2). Note that both η and r_ISCO are functions of the spacetime metric. Therefore, for given spacetime parameters, the Eddington ratio $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ can be used as the disk thickness parameter, since the geometric thickness increases as we increase $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$. In this paper, for the sake of simplicity, we will always assume that the spacetime geometry is described by the simplest version of the Johannsen metric with the only possible nonvanishing deformation parameter α₁₃, while all other deformation parameters will be set to zero (Johannsen 2013); the expression of such a metric with some basic properties is reported in the Appendix. Figure 1 illustrates our disks of finite thickness around Johannsen black holes with different values of spin and deformation parameters. We note that the pressure scale height H in Equation (1) is derived from a Newtonian model (Shakura & Sunyaev 1973). The actual pressure scale height must thus deviate from ours as we approach the black hole, affecting both the surface shape and the thickness of the disk. We can expect that the actual thickness of the disk is smaller for the same $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ and thus our model may overestimate a bit the effect; the different surface shape affects the actual value of the emission angle, and in turn, the Doppler boosting, which has an impact on the shape of the reflection spectrum. However, here the spirit is to extend our disk model from infinitesimally thin to finite thickness, as well as to get a crude estimate of the impact of the disk thickness on the measurement of the model parameters. A relativistic disk model would predict a somewhat different pressure scale height H, but to keep the analysis simple and facilitate the comparison with Taylor & Reynolds (2018a) we use their model.

Zoom In Zoom Out Reset image size

Since the spacetime is stationary and axisymmetric there are two Killing vectors, namely, a timelike and an azimuthal. Therefore, there are two conserved quantities: the specific energy E and the z-component of the specific angular momentum L_z. The system is fully determined by imposing that the gas follows nearly geodesic equatorial circular orbits (Bardeen et al. 1972).

By definition, we can write

$$\begin{eqnarray}&&\dot{t}=-\displaystyle \frac{{{Eg}}_{\phi \phi }+{L}_{z}{g}_{t\phi }}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\dot{\phi }=\displaystyle \frac{{{Eg}}_{t\phi }+{L}_{z}{g}_{{tt}}}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 3 }$$

where the overhead dot is a derivative with respect to the affine parameter (proper time for a massive particle). Employing Equations (2) and (3) in the normalization condition for the 4-velocity of massive particles u^au_a = −1, we get

$$\begin{eqnarray}&&{g}_{{rr}}{\dot{r}}^{2}+{g}_{\theta \theta }{\dot{\theta }}^{2}={V}_{\mathrm{eff}}(r,\theta ;E,{L}_{z}),\end{eqnarray} \tag{ 4 }$$

where the effective potential is

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}=-1-\displaystyle \frac{{E}^{2}{g}_{\phi \phi }+2{{EL}}_{z}{g}_{t\phi }+{L}_{z}^{2}{g}_{{tt}}}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 5 }$$

and the 4-velocity is ${u}^{a}=(\dot{t},\dot{r},\dot{\theta },\dot{\phi })$.

The explicit expressions for the energy and the angular momentum can be obtained when we impose equatorial circular orbits. The circularity condition is equivalent to require V_eff = 0 and ∂ V_eff/∂ r = 0. If we solve for E and L_z, we find

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

$$\begin{eqnarray}&&{L}_{z}=\displaystyle \frac{{g}_{t\phi }+{g}_{\phi \phi }{\rm{\Omega }}}{\sqrt{-({g}_{{tt}}+2{g}_{t\phi }{\rm{\Omega }}+{g}_{\phi \phi }{{\rm{\Omega }}}^{2})}},\end{eqnarray} \tag{ 7 }$$

where the angular velocity of equatorial circular geodesics is

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

With Ω, we can write $\dot{t}$ from ${u}^{a}{u}_{a}=-1$ and considering that for equatorial circular orbits we have ${u}^{a}=(\dot{t},0,0,\dot{\phi })\,=(1,0,0,{\rm{\Omega }})\dot{t}$,

$$\begin{eqnarray}&&\dot{t}=\displaystyle \frac{1}{\sqrt{-({g}_{{tt}}+2{g}_{t\phi }{\rm{\Omega }}+{g}_{\phi \phi }{{\rm{\Omega }}}^{2})}}.\end{eqnarray} \tag{ 9 }$$

The ISCO radius can be calculated by substituting Equations (6) and (7) into Equation (5) and solving ${\partial }^{2}{V}_{\mathrm{eff}}/\partial {r}^{2}=0$ for r.

rellxil_nk employs the formalism of the transfer function for geometrically thin and optically thick accretion disks (Cunningham 1975; Speith et al. 1995; Dauser et al. 2010). The observed reflection spectrum is the sum of the observed specific intensities I_o(ν_o) at frequency ν_o from all parts of the disk. We can perform this sum by projecting the accretion disk onto a plane perpendicular to the line of sight of the observer, which corresponds to the observer's sky (Cunningham 1975).

The observer is located at spatial infinity $(r=+\infty )$ with inclination angle ι between the normal to the disk and the line of sight of the distant observer. We use Cartesian coordinates (α, β) on the observer's plane. In terms of the photon momentum, the celestial coordinates can be written as

$$\begin{eqnarray}&&\alpha =\mathop{\mathrm{lim}}\limits_{r\to \infty }\displaystyle \frac{-{{rp}}^{(\phi )}}{{p}^{(t)}},\,\beta =\mathop{\mathrm{lim}}\limits_{r\to \infty }\displaystyle \frac{{{rp}}^{(\theta )}}{{p}^{(t)}},\end{eqnarray} \tag{ 10 }$$

where p^(a)s are the components of the 4-momentum of the photon with respect to a locally nonrotating reference frame (Bardeen et al. 1972) and are related to p^as through a coordinate transformation (e.g. p^ϕ = p^(ϕ)/ sin ι).⁷ The celestial coordinates (α, β) and the solid angle on the observer's sky are related to each other through dα dβ = D²dΩ, where D is the distance between the black hole and the observer (Cunningham 1975).

Using Liouville's theorem (Lindquist 1966), that states I_ν/ν³ = const., we can obtain the specific intensity as seen by the observer. The observed flux of an accretion disk can then be written as

$$\begin{eqnarray}&&{F}_{{\rm{o}}}({\nu }_{{\rm{o}}})=\displaystyle \frac{1}{{D}^{2}}\int {g}^{3}{I}_{{\nu }_{{\rm{e}}}}\left({r}_{{\rm{e}}},{\theta }_{{\rm{e}}}\right)d\alpha d\beta ,\end{eqnarray} \tag{ 11 }$$

where ${I}_{{\nu }_{{\rm{e}}}}({r}_{{\rm{e}}},{\theta }_{{\rm{e}}})$ is the local specific intensity, r_e is the emission radius, θ_e is the photon emission angle in the rest frame of the gas, ν_e is the photon frequency in the rest frame of the gas, and g is the redshift factor

$$\begin{eqnarray}&&g=\displaystyle \frac{{\nu }_{{\rm{o}}}}{{\nu }_{{\rm{e}}}}=\displaystyle \frac{{({p}_{a}{u}^{a})}_{{\rm{o}}}}{{({p}_{b}{u}^{b})}_{{\rm{e}}}}.\end{eqnarray} \tag{ 12 }$$

Here p^a is the 4-momentum of a photon, and ${u}_{{\rm{o}}}^{a}$ and ${u}_{{\rm{e}}}^{a}$ are the 4-velocities of the distant observer and the particles of the gas, respectively. The photon's 4-momentum is ${p}_{a}\,=(-{E}^{\gamma },{p}_{r},{p}_{\theta },{L}_{z}^{\gamma })$ and the observer is treated as static, ${u}_{{\rm{o}}}^{a}=(1,0,0,0)$. As we mentioned in Section 2, the 4-velocity of the orbiting material in the accretion disk is ${u}_{{\rm{e}}}^{a}=\dot{t}(1,0,0,{\rm{\Omega }})$, where $\dot{t}$ is given by Equation (9) and Ω is given by Equation (8). If we plug Equation (8) and Equation (9) into Equation (12), the redshift factor becomes

$$\begin{eqnarray}&&g=\displaystyle \frac{\sqrt{-({g}_{{tt}}+2{g}_{t\phi }{\rm{\Omega }}+{g}_{\phi \phi }{{\rm{\Omega }}}^{2})}}{1-{\rm{\Omega }}b},\end{eqnarray} \tag{ 13 }$$

where $b\equiv {L}_{z}^{\gamma }/{E}^{\gamma }$, which is a constant of motion along the photon trajectory.

Since the local spectrum is not isotropic, it is necessary to calculate the emission angle. The normal of the disk's surface is given by

$$\begin{eqnarray}&&{n}^{a}=\displaystyle \frac{1}{\sqrt{{g}^{{rr}}{Z}_{,r}^{2}+{g}^{\theta \theta }{Z}_{,\theta }^{2}}}(0,{g}^{{rr}}{Z}_{,r},{g}^{\theta \theta }{Z}_{,\theta },0){| }_{Z(r,\theta )},\end{eqnarray} \tag{ 14 }$$

where $Z=Z(r,\theta )$ is the "surface function" defined as

$$\begin{eqnarray}&&Z(r,\theta )=z(\rho )-2H=r\cos \theta -\displaystyle \frac{3}{\eta }\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{\mathrm{Edd}}}\right)\left[1-\sqrt{\displaystyle \frac{{r}_{\mathrm{ISCO}}}{\rho }}\right],\end{eqnarray} \tag{ 15 }$$

and goes to zero on the surface of the disk. The gradient of Z(r, θ) gives the normal to the surface. The emission angle is thus given by

$$\begin{eqnarray}&&\cos {\theta }_{{\rm{e}}}=\displaystyle \frac{g}{\sqrt{{g}^{{rr}}{Z}_{,r}^{2}+{g}^{\theta \theta }{Z}_{,\theta }^{2}}}\left[{Z}_{,r}\dot{r}+{Z}_{,\theta }\dot{\theta }\right],\end{eqnarray} \tag{ 16 }$$

where r and θ are the coordinates at the emission point in the disk.

We can define the relative redshift factor g* at a given radius of the accretion disk as (Cunningham 1975)

$$\begin{eqnarray}&&{g}^{* }=\displaystyle \frac{g-{g}_{\min }}{{g}_{\max }-{g}_{\min }}\in [0,1],\end{eqnarray} \tag{ 17 }$$

where ${g}_{\min }\,=\,{g}_{\min }({r}_{{\rm{e}}},\iota )$ and ${g}_{\max }\,=\,{g}_{\max }({r}_{{\rm{e}}},\iota )$ represent, respectively, the minimum and maximum values of the redshift factor g for the photons emitted at r_e and detected on the distant screen with inclination angle ι.

Introducing the transfer function, we can rewrite the observed flux as

$$\begin{eqnarray}&&{F}_{{\rm{o}}}({\nu }_{{\rm{o}}})=\displaystyle \frac{1}{{D}^{2}}{\int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{\int }_{0}^{1}\displaystyle \frac{\pi {r}_{{\rm{e}}}{g}^{2}f({g}^{* },{r}_{{\rm{e}}},\iota )}{\sqrt{{g}^{* }(1-{g}^{* })}}{I}_{{\rm{e}}}({r}_{{\rm{e}}},{\theta }_{{\rm{e}}}){{dg}}^{* }{{dr}}_{{\rm{e}}},\end{eqnarray} \tag{ 18 }$$

where R_out and R_in are, respectively, the outer and inner radii of the disk. Let us note that we performed a coordinate transformation from (α, β) to (r_e, g*), which means that now we carry out the integration over the accretion disk. $f({g}^{* },{r}_{{\rm{e}}},\iota )$ is the transfer function, which is given by

$$\begin{eqnarray}&&f({g}^{* },{r}_{{\rm{e}}},\iota )=\displaystyle \frac{1}{\pi {r}_{{\rm{e}}}}g\sqrt{{g}^{* }(1-{g}^{* })}\left|\displaystyle \frac{\partial (\alpha ,\beta )}{\partial ({g}^{* },{r}_{{\rm{e}}})}\right|,\end{eqnarray} \tag{ 19 }$$

where $\left|\partial (\alpha ,\beta )/\partial ({g}^{* },{r}_{{\rm{e}}})\right|$ is the Jacobian.

As noted in Taylor & Reynolds (2018a), the inner part of the accretion disk will be obscured as $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ increases. For any unobscured part of the disk, for given values of r_e and ι, the transfer function is a closed curve parametrized by g*, except in the special cases ι = 0 and π/2. There is only one point in the disk for which ${g}^{* }=0$ and one point for which ${g}^{* }=1$. There are two curves that connect these two points, so there are two branches of the transfer function, say ${f}^{(1)}({g}^{* },{r}_{{\rm{e}}},\iota )$ and ${f}^{(2)}({g}^{* },{r}_{{\rm{e}}},\iota )$. This allows us to rewrite Equation (18) as

$$\begin{eqnarray}&&\begin{array}{l}{F}_{{\rm{o}}}({\nu }_{{\rm{o}}})=\displaystyle \frac{1}{{D}^{2}}{\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{\displaystyle \int }_{0}^{1}\displaystyle \frac{\pi {r}_{{\rm{e}}}{g}^{2}{f}^{(1)}({g}^{* },{r}_{{\rm{e}}},\iota )}{\sqrt{{g}^{* }(1-{g}^{* })}}{I}_{{\rm{e}}}({r}_{{\rm{e}}},{\theta }_{{\rm{e}}}^{(1)}){{dg}}^{* }{{dr}}_{{\rm{e}}}\\ \ +\ \displaystyle \frac{1}{{D}^{2}}{\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{\displaystyle \int }_{0}^{1}\displaystyle \frac{\pi {r}_{{\rm{e}}}{g}^{2}{f}^{(2)}({g}^{* },{r}_{{\rm{e}}},\iota )}{\sqrt{{g}^{* }(1-{g}^{* })}}{I}_{{\rm{e}}}({r}_{{\rm{e}}},{\theta }_{{\rm{e}}}^{(2)}){{dg}}^{* }{{dr}}_{{\rm{e}}},\end{array}\end{eqnarray} \tag{ 20 }$$

where ${\theta }_{{\rm{e}}}^{(1)}$ and ${\theta }_{{\rm{e}}}^{(2)}$ present the emission angles with relative redshift factor g* in the branches 1 and 2, respectively. For values of r_e for which the disk is obscured, some portion of ${f}^{(1)}({g}^{* },{r}_{{\rm{e}}},\iota )$ and/or ${f}^{(2)}({g}^{* },{r}_{{\rm{e}}},\iota )$ corresponding to the obscured parts of the disk will be equal to zero (see Figures 2–4), which means there is no radiation contributing from this part into the total reflection spectrum. In such cases, the integration in Equation (20) is performed only on nonzero values of the transfer function.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Numerical Method

Here we describe our method for calculating the transfer function and, following the methodology of relxill and relxill_nk, we create a FITS (Flexible Image Transport System) file that contains the relevant spacetime information. The structure of the FITS file is similar to that used in relxill_nk for a infinitesimally thin accretion disk and is fully described in Abdikamalov et al. (2019). There are three physical parameters describing the black hole spacetime in the table, namely, the dimensionless black hole spin parameter a_*, the deformation parameter (α₁₃ for the metric considered in this paper), and the inclination angle ι, structured in a 30 by 30 by 22 grid, respectively. The grid for the black hole spin becomes denser as the black hole spin increases, since for high values of a_* the ISCO changes faster. The values of the deformation parameters of the Johannsen metric in the grid are first evenly distributed in the range [−5, 5]. However, for negative values of α₁₃ we may have spacetimes with pathological properties; see the Appendix and the constraint on α₁₃ in Equation (A4). The lower bound of α₁₃ is thus set to the larger value between −5 and the bound in Equation (A4). The grid point along the inclination angles are evenly distributed in 0 < cos ι < 1. For each combination of a_*, α₁₃, and ι, we discretize the accretion disk into 100 emission radii r_e. For every emission radius, the transfer functions, f, and emission angles, θ_e, are tabulated at 40 equally spaced values of g* on each branch of the transfer function.

As in Abdikamalov et al. (2019), a general relativistic ray-tracing code is used to compute the necessary parameters for the FITS file, namely, the redshift factor, emission angle, and the Jacobian. The ray-tracing code calculates the trajectories of photons in the Johannsen metric from the black hole accretion disk to a distant observer. The code follows the method described in Psaltis & Johannsen (2012) and is a modified version of the one used in Abdikamalov et al. (2019), Ayzenberg & Yunes (2018), Gott et al. (2019). Since it is a stationary and axisymmetric spacetime, the Johannsen metric has a conserved energy E and a conserved angular momentum L_z. Their relation to the 4-momentum of a test particle, p_t = −E and p_ϕ = L_z, leads to two first-order differential equations shown in Equations (2) and (3). Rewriting these two equations in terms of the impact parameter b ≡ L_z/E and the normalized affine parameter $\lambda ^{\prime} \equiv E\lambda$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{dt}}{d\lambda ^{\prime} }=-\displaystyle \frac{{{bg}}_{t\phi }+{g}_{\phi \phi }}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{d\lambda ^{\prime} }=\displaystyle \frac{{g}_{t\phi }+{{bg}}_{{tt}}}{{g}_{{tt}}{g}_{\phi \phi }-{g}_{t\phi }^{2}}.\end{eqnarray} \tag{ 22 }$$

The second-order geodesic equations for a generic axisymmetric metric describe the evolution of the r- and θ-components of the photon's position as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}r}{d\lambda {{\prime} }^{2}}=-{{\rm{\Gamma }}}_{{tt}}^{r}{\left(\displaystyle \frac{{dt}}{d\lambda ^{\prime} }\right)}^{2}-{{\rm{\Gamma }}}_{{rr}}^{r}{\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)}^{2}-{{\rm{\Gamma }}}_{\theta \theta }^{r}{\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right)}^{2}\\ -{{\rm{\Gamma }}}_{\phi \phi }^{r}{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}^{2}-2{{\rm{\Gamma }}}_{t\phi }^{r}\left(\displaystyle \frac{{dt}}{d\lambda ^{\prime} }\right)\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)\\ -2{{\rm{\Gamma }}}_{r\theta }^{r}\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right),\end{array}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}\theta }{d\lambda {{\prime} }^{2}}=-{{\rm{\Gamma }}}_{{tt}}^{\theta }{\left(\displaystyle \frac{{dt}}{d\lambda ^{\prime} }\right)}^{2}-{{\rm{\Gamma }}}_{{rr}}^{\theta }{\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)}^{2}-{{\rm{\Gamma }}}_{\theta \theta }^{\theta }{\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right)}^{2}\\ -{{\rm{\Gamma }}}_{\phi \phi }^{\theta }{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}^{2}-2{{\rm{\Gamma }}}_{t\phi }^{\theta }\left(\displaystyle \frac{{dt}}{d\lambda ^{\prime} }\right)\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)\\ -2{{\rm{\Gamma }}}_{r\theta }^{\theta }\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right),\end{array}\end{eqnarray} \tag{ 24 }$$

where ${{\rm{\Gamma }}}_{{bc}}^{a}$ indicate the Christoffel symbols of the metric.

A coordinate system and reference frame are chosen in such a way that the black hole is located at the origin and the black hole spin angular momentum is along the z-axis. We set M = 1 in what follows and in the code, since the reflection spectrum does not directly depend on the black hole mass M. The observer screen is located far from the black hole at a distance of D = 10⁸, with azimuthal angle θ = ι and polar angle ϕ = 0. The screen uses the polar coordinates r_scr and ϕ_scr, and their relation to the celestial coordinates of Equation (10) are $\alpha ={r}_{\mathrm{scr}}\cos {\phi }_{\mathrm{scr}}$ and $\beta ={r}_{\mathrm{scr}}\sin {\phi }_{\mathrm{scr}}$.

The code solves the system of equations—Equations (21)–(24)—backwards in time. Each photon has an initial position on the screen and an initial 4-momentum perpendicular to the screen. In the Boyer–Lindquist coordinates of the black hole spacetime, the initial position and 4-momentum of the photon are given by

$$\begin{eqnarray}&&{r}_{{\rm{i}}}={\left({\alpha }^{2}+{\beta }^{2}+{D}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\theta }_{{\rm{i}}}=\arccos \left(\displaystyle \frac{D\cos \iota +\beta \sin \iota }{{r}_{{\rm{i}}}}\right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\phi }_{{\rm{i}}}=\arctan \left(\displaystyle \frac{\alpha }{D\sin \iota -\beta \cos \iota }\right),\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)}_{{\rm{i}}}=\displaystyle \frac{D}{{r}_{{\rm{i}}}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}=\displaystyle \frac{-\cos \iota +\tfrac{D}{{r}_{{\rm{i}}}^{2}}\left(D\cos \iota +\beta \sin \iota \right)}{\sqrt{{r}_{{\rm{i}}}^{2}-{\left(D\cos \iota +\beta \sin \iota \right)}^{2}}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}=\displaystyle \frac{-\alpha \sin \iota }{{\alpha }^{2}+{\left(D\sin \iota -\beta \cos \iota \right)}^{2}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{{dt}}{d\lambda ^{\prime} }\right)}_{{\rm{i}}}=-\displaystyle \frac{{g}_{t\phi }}{{g}_{{tt}}}{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}\\ \ +\ \sqrt{\displaystyle \frac{{g}_{t\phi }^{2}}{{g}_{{tt}}^{2}}{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}^{2}-\left[\displaystyle \frac{{g}_{{rr}}}{{g}_{{tt}}}{\left(\displaystyle \frac{{dr}}{d\lambda ^{\prime} }\right)}_{{\rm{i}}}^{2}+\displaystyle \frac{{g}_{\theta \theta }}{{g}_{{tt}}}{\left(\displaystyle \frac{d\theta }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}^{2}+\displaystyle \frac{{g}_{\phi \phi }}{{g}_{{tt}}}{\left(\displaystyle \frac{d\phi }{d\lambda ^{\prime} }\right)}_{{\rm{i}}}^{2}\right]}.\end{array}\end{eqnarray} \tag{ 31 }$$

The expression for ${\left({dt}/d\lambda ^{\prime} \right)}_{{\rm{i}}}$ is obtained by requiring that the norm of the photon 4-momentum vanishes. The impact parameter b is a conserved quantity used in Equations (21) and (22), calculated from the initial conditions.

Our algorithm adaptively searches for the photons that hit the surface of the accretion disk, i.e. the z(r_e) surface, at the 100 disk emission radii r_e to within a precision of ∼10⁻⁷ varying r_scr and ϕ_scr. For each emission radius we first shoot 10 photons, from which we register preliminary g_min and g_max. Then the actual redshift extremas are found by shifting ϕ_scr from these preliminary extremas with an adaptive step-size. Afterwards we search for 80 different photons, 40 in each branch of the transfer function, that correspond to equally distributed values of g*∈[0,1]. The photons are split into two branches according to

$$\begin{eqnarray}&&{\phi }_{\mathrm{scr}}^{\min }\lt {\phi }_{\mathrm{scr}}\lt {\phi }_{\mathrm{scr}}^{\max }\ \mathrm{and}\ {\phi }_{\mathrm{scr}}^{\min }\gt {\phi }_{\mathrm{scr}}\gt {\phi }_{\mathrm{scr}}^{\max }\end{eqnarray} \tag{ 32 }$$

where ${\phi }_{\mathrm{scr}}^{\min }$ and ${\phi }_{\mathrm{scr}}^{\max }$ are the photons with g_min and g_max, respectively. When searching for photons at a disk emission radius ${r}_{{\rm{e}}}$, we divide photons into real and imaginary ones. Imaginary photons are those that cross the disk several times before landing on the target ring of the disk. Similarly, a real photon does not cross the disk before landing on the target radius. This separation helps to distinguish photons originating from obscured and unobscured parts of the disk. Therefore, imaginary photons originate from the obscured part of the disk that we cannot see.

For each of these photons, we calculate the redshift factor g, Equation (12), emission angle θ_e, Equation (16), and Jacobian $| \partial (\alpha ,\beta )/\partial ({g}^{* },{r}_{{\rm{e}}})|$. The latter is calculated by using

$$\begin{eqnarray}&&\left|\displaystyle \frac{\partial (\alpha ,\beta )}{\partial ({g}^{* },{r}_{{\rm{e}}})}\right|={\left(-1\right)}^{p}\left({g}_{\max }-{g}_{\min }\right)\left[\displaystyle \frac{\partial \alpha }{\partial g}\displaystyle \frac{\partial \beta }{\partial {r}_{{\rm{e}}}}-\displaystyle \frac{\partial \alpha }{\partial {r}_{{\rm{e}}}}\displaystyle \frac{\partial \beta }{\partial g}\right],\end{eqnarray} \tag{ 33 }$$

where p is the number of the branch (1 for branch 1 and 2 for branch 2). The introduction of imaginary photons avoids the zero values of the Jacobian, hence, the zero values of the transfer function. We attribute a negative sign to them because the negative transfer function allows us to leave unchanged the standard interpolation scheme used in relxill_nk. We set negative values of the transfer function to zero only after the final interpolation in the model. The third term on the right-hand side is computed using an adaptive algorithm that, when solving the geodesic equations, searches for two photons that have g ± Δg for the given emission radius r_e, where g is the initial redshift factor of the original and Δg < 10⁻⁶. For r_e ± Δr, the code uses adjacent photons from two neighboring emission rings, thus, Δr is the distance between these photons. The derivatives are then calculated from the emission radius, redshift factor, and initial coordinates of these four photons in a separate code, as the second term on the right-hand side of Equation (33).

Finally, we use a separate script to process all photons and create the FITS file for a specific value of $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$. First, the script calculates the Jacobian and then generates a FITS file containing the values of 100 emission radii r_e, corresponding minimum and maximum redshift values (g_min and g_max), transfer functions, and emission angles θ_e, for the full combination of dimensionless spin a_*, deformation parameter α₁₃, and inclination angle ι.

3.2. Impact of Disk Thickness on the Iron Lines

Figures 5–7 show iron line profiles for different values of a_*, α₁₃, ι, and $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$. All calculations assume that the energy of the line in the rest frame of the gas is E = 6.4 keV and that the intensity profile of the disk is described by a power law with emissivity index q = 3, i.e., I_e  ∝ ρ⁻³. Here the intensity profile is described using the coordinate ρ, and not r, because such a choice guarantees a one-to-one correspondence between every value of ρ and every point of the surface of the disk (modulo the cylindrical symmetry of the system). At a qualitative level, we can say that the impact of the thickness of the disk on the iron line profile is weak for ι = 20° and 45° (Figures 5 and 6), while it is a bit larger for ι = 70° (Figure 7). For a_* and α₁₃ there is not a clear trend: generally speaking, if the ISCO radius is closer to the black hole, the gravitational field is stronger, and we can expect that small variations in the exact emission point has a larger impact on the shape of the iron line profile; however, when the ISCO radius is closer to the black hole, η is typically higher, making the disk thinner for the same value of $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$, thus producing the opposite effect with smaller difference with respect to an infinitesimally thin disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 7, we see a peaky feature in some iron line profiles for disks of finite thickness, while such a feature is never present in the case of infinitesimally thin disks. It is not a numerical artifact but the result of a combination of the disk self-shadowing, namely, the fact that a disk of finite thickness can obscure a portion of the very inner part of the disk and relativistic effects. We note that the effect of self-shadowing shows up above some critical value of the viewing angle (which depends on the black hole spin, deformation parameter, mass accretion rate, and emissivity profile) and the effect is more and more important as the viewing angle increases, as both the self-shadowing and the Doppler boosting increase. In the Kerr spacetime (α₁₃ = 0) we see a peaky feature in Figure 7, where ι = 70°, and for a_* = 0.9. When a_* = 0.998, the peaky feature disappears because the radiative efficiency η is significantly higher and the disk is thus thinner, see Equation (1) and Figure 1. Taylor & Reynolds (2018a) do not find such a peaky feature in their paper because their plots show iron line profiles observed with viewing angles up to 60°, while Figure 7 is for 70°. Moreover, the feature depends also on the disk intensity profile: here we use a power-law profile, while Taylor & Reynolds (2018a) employ a lamppost emissivity profile.

To illustrate the origin of the peaky feature, in Figure 8 we show the iron line profile contributions from a few inner annuli of the accretion disk for $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0$ (left panel) and $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.3$ (right panel) in the case a_* = 0.9, ι = 70°, and q = 3. As can be seen from the comparison of the two panels, the spectra of the inner annuli (r_ISCO < ρ < 3 M) is significantly different in the two iron line profiles, while the spectra of outer annuli are quite similar. The spectra of annuli with r > 10 M are not shown in Figure 8 because the difference is negligible. While the spectrum of the inner annulus of the infinitesimally thin disk has a broad peak, that of the disk with finite thickness has a sharp cutoff at high energy due to the obscuration of the disk. The Doppler effect and light bending are also different as the gas keeps the same angular velocity but leaves the equatorial plane of the infinitesimally thin disk; indeed, the emission angle θ_e changes and the inner part of the accretion disk behind the black hole becomes more face on as the thickness of the disk increases, affecting the Doppler boosting.

Zoom In Zoom Out Reset image size

Fast-rotating black holes observed from large viewing angles are the most interesting and suitable sources for testing the Kerr metric, since these two properties maximize the relativistic effects in the reflection spectrum. As we have seen in this section, these turn out to be even the sources with the reflection spectra more affected by the thickness of the disk. It is thus important to further investigate the features of their spectra. In particular, a relevant parameter is the emissivity profile. In Figures 5–7, we have calculated iron line profiles employing an emissivity index q = 3. In the analysis of real data, it is common to find much higher emissivity indices, especially for the inner part of the disk (an example is presented in the next section). If we increase the value of q, we increase the fraction of radiation emitted from the very inner part of the accretion disk, enhancing the impact of the thickness of the disk on the iron line shape.

Figure 9 shows iron line shapes for a spin parameter a_* = 0.9, a viewing angle ι = 70° and 80°, and an emissivity index q = 10 for infinitesimally thin disks and disks of finite thickness. Figure 10 is as Figure 9 but for a_* = 0.998. We note that employing q = 6 does not change the line shapes much, while there are clear differences for lower values of q. While the effects of a_*, α₁₃, ι, q, and $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ mix together in a quite complicated way, we can see that the main impact of the disk thickness is to move the peak of the iron line shape to lower energies. This is again the effect of the obscuration of the inner part of the accretion disk and of the different Doppler boosting discussed a few paragraphs above and related to the peaky feature of the disk with finite thickness. Enhancing the fraction of radiation emitted from the very inner part of the accretion disk by increasing the value of q leads to increase the effect of the thickness of the disk on the shape of the iron line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the FITS file for a specific value of $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$, we can test the new model with real data in order to estimate the systematic uncertainties of the model with an infinitesimally thin disk and the impact of the disk thickness on our tests of the Kerr black hole hypothesis. For a preliminary study to present in this paper, we consider the 2007 Suzaku observation of GRS 1915+105, which was analyzed for the first time by Blum et al. (2009) with a Kerr model and was analyzed by our group in Zhang et al. (2019a) to test the Kerr metric with relxill_nk. These data provide, as of now, one of the most stringent constraints on α₁₃ among all the observations and sources analyzed so far with relxill_nk, and they are thus suitable to test the impact of the disk thickness on the measurement of the deformation parameter α₁₃. Note that during the 2007 Suzaku observation, the Eddington-scaled accretion luminosity of the source was about 0.2,⁸ which is a high value but still in the range expected for a geometrically thin accretion disk with inner edge at the ISCO radius.

The observation, data reduction, and choice of the model has already been discussed in Zhang et al. (2019a). Here we just point out the main properties of this observation and this source. GRS 1915+105 is quite a bright stellar-mass black hole. Previous analyses in Kerr and Johannsen backgrounds suggest that the inclination angle of the disk is high (around 70°) and the inner edge of the accretion disk is very close to the black hole. A high value of the inclination angle tends to maximize the relativistic effects, in particular the light bending. Moreover, as we have seen at the end of the previous section, it seems that a high inclination angle maximizes the impact of the thickness of the accretion disk, which is indeed what we want to explore here. Concerning the position of the inner edge of the accretion disk, we meet two opposing effects. If the inner edge is closer to the black hole, the signature of the strong gravity effects in the reflection spectrum are larger, and it may be possible that small differences in the location of the emission can have an impact on the reflection spectrum. Note also that in our tests of the Kerr black hole hypothesis we typically prefer to analyze sources with an inner edge of the accretion disk very close to the compact object because this helps to break the parameter degeneracy and constrain the deformation parameter; if the inner edge of the accretion disk is far from the source, simultaneous measurements of the black hole spin and the deformation parameter are difficult or impossible. On the contrary, an inner edge of the accretion disk very close to the black hole is typically accompanied by a high value of the radiative efficiency η, which makes the disk thinner via Equation (1). The difference with the infinitesimally thin disk should thus be smaller.

The analysis of the 2007 Suzaku observation of GRS 1915+105 reported in Zhang et al. (2019a) provides, at the moment, one of the most stringent and robust tests of the Kerr metric with relxill_nk. Suzaku observed GRS 1915+105 on 27 May 2007 (obs. ID 402071010), when the source was in the low-hard state, for approximately 117 ks. After all efficiencies and screening, the net exposure time is 29 ks for the XIS1 camera and 53 ks for HXD/PIN, while the other XIS units either were turned off or run in a special timing mode. As shown in Zhang et al. (2019a), the hardness of the source was quite stable in the 2007 Suzaku observation. The spectrum is clearly dominated by a strong relativistic reflection component, with a clear broad iron line around 6 keV and a Compton hump peaking around 20 keV. We do not see any thermal component from the disk, which is also welcome because the nonrelativistic reflection model employed is xillver, which should only be used for cold accretion disks. The quality of the Suzaku data is very good and we have both a high energy resolution near the iron line with the XIS1 instrument and a broad energy band when we add the PIN data.

As discussed in Zhang et al. (2019a), the XSPEC model tbabs×relxill_nk fits the data well and it seems that we do not need other components.⁹ tbabs describes the Galactic absorption (Wilms et al. 2000), and the hydrogen column density is left free in all fits. relxill_nk is our relativistic reflection spectrum in the Johannsen metric with nonvanishing deformation parameter α₁₃, and we consider two models: the infinitesimally thin disk ($\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0$) and the model with $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$. The best-fit values for the two models are reported in Table 1, where the parameter uncertainties correspond to the 90% confidence level. Best-fit models and ratio plots are shown in Figure 11. As our interest here is in the impact of the disk thickness on tests of the Kerr metric, we also show the constraints on the black hole spin and the deformation parameter in Figure 12 after marginalization over all other free parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Best-fit Values from the Analysis of the 2007 Suzaku Sbservation of GRS 1915+105 relxill_nk Employing an Infinitesimally Thin Disk (Left Column) and a Disk with Finite Thickness for $\dot{M}/{\dot{M}}_{{\rm{E}}{\rm{d}}{\rm{d}}}=0.2$

	Infinitesimally Thin Disk	Disk with Finite Thickness
tbabs
${N}_{{\rm{H}}}/{10}^{22}$ cm−2	${7.97}_{-0.09}^{+0.07}$	${7.867}_{-0.024}^{+0.022}$
relxill_nk
qin	10.0−0.6	${8.55}_{-1.01}^{+0.13}$
qout	0.00+0.21	0.0+1.1
Rbr [M]	${6.03}_{-0.44}^{+0.18}$	${7.26}_{-0.11}^{+3.62}$
i [deg]	${73.7}_{-0.6}^{+1.6}$	${79.6}_{-0.5}^{+3.3}$
a*	${0.9897}_{-0.0009}^{+0.0015}$	${0.9950}_{-0.0003}^{({\rm{P}})}$
α13	$-{0.09}_{-0.10}^{+0.10}$	${0.00}_{-0.15}^{+0.05}$
$\dot{M}/{\dot{M}}_{\mathrm{Edd}}$	0a	0.2a
$\mathrm{log}\xi$	${2.77}_{-0.04}^{+0.03}$	${2.699}_{-0.010}^{+0.011}$
${A}_{\mathrm{Fe}}$	${0.60}_{-0.06}^{+0.07}$	${0.737}_{-0.032}^{+0.021}$
Γ	${2.199}_{-0.016}^{+0.015}$	${2.2120}_{-0.0016}^{+0.0059}$
${E}_{\mathrm{cut}}$ [keV]	${71.2}_{-1.6}^{+3.3}$	${69.6}_{-1.1}^{+0.5}$
${R}_{{\rm{f}}}$	${0.48}_{-0.03}^{+0.09}$	${0.461}_{-0.073}^{+0.006}$
norm	${0.0429}_{-0.0025}^{+0.0004}$	${0.04626}_{-0.0044}^{+0.0005}$
${\chi }^{2}/\nu$	2314.75/2208	2306.52/2208
	=1.04835	=1.04462

Note. The reported uncertainties correspond to the 90% confidence level for one relevant parameter (${\rm{\Delta }}{\chi }^{2}=2.71$).

^aIndicates that the parameter is frozen in the fit. Note that ${q}_{\mathrm{in}}$ and ${q}_{\mathrm{out}}$ are allowed to vary in the range [0, 10] and the best fits are stuck at the boundary with the exception of ${q}_{\mathrm{in}}$ for the model with $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$. The maximum value of the spin parameter allowed by the model is 0.998, and for $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$ the 90% confidence level reaches the boundary.

Download table as:  ASCIITypeset image

The fit of the model with a disk of finite thickness is only a bit better, but not significantly better than the model with an infinitesimally thin disk (${\rm{\Delta }}{\chi }^{2}=8.23$). The measurements of most model parameters are consistent; in particular, there is no difference in the final constraint on the deformation parameter α₁₃. Both models require a very high spin parameter; the two measurements are slightly different if we believe in the statistical uncertainty of the fits, but that is indeed too low to expect that systematic uncertainties are not dominant. The model parameter presenting some difference in the two measurements is the disk inclination angle, and the model with a disk of finite thickness requires a very high value of ι.

In both models, we find a very high inner emissivity index and a very low outer emissivity index. We interpret this result as a possible indication of a corona with a ring-like axisymmetry geometry located just above the accretion disk, which actually would be fitted better with a twice broken power law with very steep emissivity profile over the inner region, then flattening in the intermediate region, and falling off approximately as r⁻³ over the outer region (Miniutti et al. 2003; Wilkins & Fabian 2011; Wilkins & Gallo 2015b). Such a coronal geometry above the accretion disk would predict the Comptonization of the relativistic reflection component that, when not taken into account in the XSPEC model, would lead to residuals similar to those shown in the lower panels of Figure 11 (Wilkins & Gallo 2015a). We note that other authors have interpreted very high inner emissivity indices from the fit as a deficiency of the model, in particular of the assumption of a constant ionization profile of the disk (Kammoun et al. 2019). Such an interpretation would presumably lead to a different measurement of the spin and of the deformation parameter (Shreeram & Ingram 2020), and we plan to leave the study of such a possibility and its impact in the constraint on α₁₃ to future work.

Since the model tbabs×relxill_nk has 13 free parameters in our fits and the χ² minimizing algorithm of XSPEC has often problems to reliably find a minimum in complicated χ² landscapes, we perform a Markov Chain Monte Carlo (MCMC) analysis of the case with $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$ using the Python script by Jeremy Sanders, which uses emcee (MCMC Ensemble sampler implementing the Goodman & Weare algorithm).¹⁰ We use 200 walkers of 10,000 iterations, after burning the initial 1000 iterations (which is around 100 times the autocorrelation length). Figure 13 shows the corner plot with all the one- and two-dimensional projections of the posterior probability distributions of all the free parameters (we only omit the normalization of relxill_nk). The two-dimensional projections report the 1, 2, and 3σ confidence level limits for two relevant parameters. In Figure 14 we zoom into the spin parameter versus deformation parameter panel of Figure 13. We note that the results of the MCMC analysis is consistent with the results obtained with XSPEC; in particular, the constraints in Figure 14 are very similar to the constraints in the right panel in Figure 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The possibility of performing precision tests of general relativity in the strong gravity region around black holes using X-ray reflection spectroscopy is determined by our capability of limiting the systematic uncertainties (broadly defined) in the final measurement of possible deviations from the Kerr background. The work presented in this paper is a step of our program to develop a sufficiently sophisticated relativistic reflection model to perform precision tests of the Kerr black hole hypothesis.

In current relativistic reflection models, the accretion disk is assumed to be infinitesimally thin, while in reality it has a finite thickness, which should increase as the mass accretion rate increases. Here we have presented an extension of relxill_nk in which the disk has a finite thickness by implementing the disk geometry proposed in Taylor & Reynolds (2018a). With the current structure of the model, we cannot add the mass accretion rate as a new model parameter capable of varying over some range, as this would make the FITS file too large. We have thus constructed FITS files for specific values of the mass accretion rate of the source. With our current version of the ray-tracing code, the construction of a single FITS file for a specific value of $\dot{M}/{\dot{M}}_{\mathrm{Edd}}$ requires about 2 weeks on a computer cluster with about 250 cores. The size of the FITS file is about 1.3 GB.

In Section 4 we have analyzed the 2007 Suzaku observation of GRS 1915+105 with relxill_nk assuming either that the accretion disk is infinitesimally thin and that the disk has a finite thickness with $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$, which is the estimate inferred from the Suzaku observation and the known values of mass and distance of the source. Our analysis does not show a significant difference in the estimate of the model parameters, and in particular, in the constraint on the deformation parameter α₁₃. We should stress that we have analyzed very high-quality data: GRS 1915+105 is a bright source and Suzaku has both a good energy resolution near the iron line (which is the most informative part of the reflection spectrum concerning the spacetime metric) and high energy data to fit the Compton hump. The source is also characterized by a high disk inclination angle, which should maximize the impact of the thickness of the disk. As of now, the analysis of these data provides one of the most stringent constraints on the Kerr metric with relxill_nk, so this motivated us to use the new model with this observation. It is possible that the weak impact of the disk thickness on the analysis of this source is related to the fact that the estimate of the radiative efficiency η is high, which makes the disk quite thin even if $\dot{M}/{\dot{M}}_{\mathrm{Edd}}=0.2$. However, this is always the case for sources used to test the Kerr metric, because for low values of η the signature of the background metric on the reflection spectrum is weak and we cannot constrain the deformation parameter due to parameter degeneracy.

Last, we note that there is no disagreement between our results and those found in Taylor & Reynolds (2018a), but a comparison is not straightforward. In our model the intensity profile is described a broken power law, while Taylor & Reynolds (2018a) consider the profile generated by a point-like lamppost corona. Taylor & Reynolds (2018a) find that the disk thickness leads to underestimating the black hole spin parameter when the data are fitted with a infinitesimally thin disk model, but their input parameter is a_* = 0.9, so η is lower and the thickness of the disk is higher. Moreover, they assume a point-like lamppost corona with height h = 3 M: for such a low value of h, the difference of the intensity profile between a disk of finite thickness and an infinitesimally thin disk is quite pronounced. In our case, since we have analyzed a source with high η, the thickness of the disk is lower and probably for this reason we do not see any clear modeling bias in the measurements of the model parameters. Note that the purpose of implementing a disk with finite thickness in relxill_nk is not primarily to fit sources with thicker disk. Our goal is to get stringent constraints on the deformation parameters and for this reason we have analyzed the Suzaku observation of GRS 1915+105, as it represents one of the most stringent tests of the Kerr metric.

We want to stress that the thickness of the disk is one of the current common model simplifications among many others, and presumably not the most crucial one. The present extension of relxill_nk implementing a disk with finite thickness does not aim to perform generic tests of the Kerr metric using X-ray reflection spectroscopy. It is meant to arrive at an estimate of the impact of the disk thickness on our capabilities of testing the Kerr black hole hypothesis. The model has still a number of simplifications that inevitably lead to modeling bias currently not well under control. Uncertainties in the coronal geometry, simplifications in the calculations of the reflection spectrum in the rest frame of the accreting gas on the disk, impact of magnetic fields on the disk structure, etc. are all effects that need to be investigated in order to improve the capability of X-ray reflection spectroscopy to study accreting black holes.

This work was supported by the Innovation Program of the Shanghai Municipal Education Commission, Grant No. 2019-01-07-00-07-E00035, and the National Natural Science Foundation of China (NSFC), Grant No. 11973019. A.B.A. also acknowledges support from the Shanghai Government Scholarship (SGS). J.A.G. and S.N. acknowledge support from the Alexander von Humboldt Foundation. C.B., J.A.G., S.N., and A.T. are members of the International Team 458 at the International Space Science Institute (ISSI), Bern, Switzerland, and acknowledge support from ISSI during the meetings in Bern.

In Boyer–Lindquist-like coordinates, the line element of the Johannsen metric reads (Johannsen 2013)

$$\begin{eqnarray}&&\begin{array}{l}{{ds}}^{2}=-\displaystyle \frac{\tilde{{\rm{\Sigma }}}\left({\rm{\Delta }}-{a}^{2}{A}_{2}^{2}{\sin }^{2}\theta \right)}{{B}^{2}}{{dt}}^{2}+\displaystyle \frac{\tilde{{\rm{\Sigma }}}}{{\rm{\Delta }}}{{dr}}^{2}+\tilde{{\rm{\Sigma }}}d{\theta }^{2}\\ -\displaystyle \frac{2a\left[\left({r}^{2}+{a}^{2}\right){A}_{1}{A}_{2}-{\rm{\Delta }}\right]\tilde{{\rm{\Sigma }}}{\sin }^{2}\theta }{{B}^{2}}{dtd}\phi \\ \ +\ \displaystyle \frac{\left[{\left({r}^{2}+{a}^{2}\right)}^{2}{A}_{1}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta \right]\tilde{{\rm{\Sigma }}}{\sin }^{2}\theta }{{B}^{2}}d{\phi }^{2},\end{array}\end{eqnarray} \tag{ A1 }$$

where M is the black hole mass, $a=J/M$, J is the black hole spin angular momentum, $\tilde{{\rm{\Sigma }}}={\rm{\Sigma }}=f$, and

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }} & = & {r}^{2}+{a}^{2}{\cos }^{2}\theta ,\\ {\rm{\Delta }} & = & {r}^{2}-2{Mr}+{a}^{2},\\ B & = & \left({r}^{2}+{a}^{2}\right){A}_{1}-{a}^{2}{A}_{2}{\sin }^{2}\theta .\end{array}\end{eqnarray} \tag{ A2 }$$

The functions f, A₁, A₂, and A₅ are defined as

$$\begin{eqnarray}\begin{array}{rcl}f & = & \displaystyle \sum _{n=3}^{\infty }{\epsilon }_{n}\displaystyle \frac{{M}^{n}}{{r}^{n-2}},\\ {A}_{1} & = & 1+\displaystyle \sum _{n=3}^{\infty }{\alpha }_{1n}{\left(\displaystyle \frac{M}{r}\right)}^{n},\\ {A}_{2} & = & 1+\displaystyle \sum _{n=2}^{\infty }{\alpha }_{2n}{\left(\displaystyle \frac{M}{r}\right)}^{n},\\ {A}_{5} & = & 1+\displaystyle \sum _{n=2}^{\infty }{\alpha }_{5n}{\left(\displaystyle \frac{M}{r}\right)}^{n},\end{array}\end{eqnarray} \tag{ A3 }$$

where {_n}, {α_1n}, {α_2n}, and {α_5n} are four infinite sets of deformation parameters without constraints from the Newtonian limit and solar system experiments. The leading order deformation parameters are thus ₃, α₁₃, α₂₂, and α₅₂. In this paper, we have only considered the deformation parameter α₁₃ because it has the strongest impact on the shape of the reflection spectrum, but all our results can be easily extended to metrics with other nonvanishing deformation parameters as well as, more in general, to any stationary, axisymmetric, and asymptotically flat black hole spacetime.

In order to avoid spacetimes with pathological properties, we must impose some constraints on the values of a_* and α₁₃. As in the case of the Kerr spacetime, we must impose that $| {a}_{* }| \leqslant 1;$ for $| {a}_{* }| \gt 1$ there is no event horizon and the solution describes the spacetime of a naked singularity. As discussed in Tripathi et al. (2018), we have to impose the following constraint on α₁₃:

$$\begin{eqnarray}&&{\alpha }_{13}\gt -\displaystyle \frac{1}{2}{\left(1+\sqrt{1-{a}_{* }^{2}}\right)}^{4}.\end{eqnarray} \tag{ A4 }$$