Reflection spectra of thick accretion disks

TL;DR

This paper investigates how thick accretion disks around black holes affect X-ray reflection spectra, challenging the thin disk assumption and exploring implications for black hole parameter estimation.

Contribution

It introduces the Polish donut model for thick disks and analyzes how disk geometry influences reflection spectra and parameter inference.

Findings

Thick disks produce distinct iron line shapes.

Disk thickness impacts the accuracy of Kerr hypothesis tests.

Simulated spectra reveal biases when using thin disk models for thick disks.

Abstract

Relativistic reflection features are commonly observed in the X-ray spectra of stellar-mass and supermassive black holes and originate from illumination of the inner part of the accretion disk by a hot corona. All the available relativistic reflection models assume that the disk is infinitesimally thin and the inner edge is at the innermost stable circular orbit or at a larger radius. However, we know that several sources, especially among supermassive black holes, have quite high mass accretion rates. In such a case, the accretion disk becomes geometrically thick and the inner edge of the disk is expected to be inside the innermost stable circular orbit. In this work, we employ the Polish donut model to describe geometrically thick disks and we study the iron line shapes from similar systems. We also simulate full reflection spectra and we analyze the simulated observations with a thin disk relativistic reflection model to determine the impact of the disk structure on the estimation of the model parameters, in particular in the case of tests of the Kerr hypothesis.

Figures

Tables

Table 1. Table 1: Input parameters and best-fit values for simulations A-D. The reported uncertainties correspond to 90% confidence level for one relevant parameter. ⋆ indicates that the parameter is frozen in the fit.

	Simulation A			Simulation B			Simulation C			Simulation D
	Input	Fit		Input	Fit		Input	Fit		Input	Fit
tbabs
$N_H/{10}^{20}$ cm-2	$6.74$	${6.74}^{\star }$		$6.74$	${6.74}^{\star }$		$6.74$	${6.74}^{\star }$		$6.74$	${6.74}^{\star }$
relxill_nk
$q$	$9$	${4.979}_{-0.010}^{+0.011}$		$9$	${5.306}_{-0.034}^{+0.011}$		$9$	${5.34}_{-0.05}^{+0.04}$		$9$	${8.28}_{-0.21}^{+0.09}$
$i$ [deg]	$35$	${27.05}_{-0.18}^{+0.14}$		$35$	${31.38}_{-0.14}^{+0.24}$		$35$	${31.7}_{-0.3}^{+0.5}$		$35$	${46.3}_{-0.4}^{+3.5}$
$a_{\ast }$	$0.80$	${0.9251}_{-0.0011}^{+0.0010}$		$0.85$	${0.9727}_{-0.0012}^{+0.0008}$		$0.90$	${0.9910}_{-0.0011}^{+0.0009}$		$0.95$	$>0.997$
${\alpha }_{13}$	$0$	$-{1.08}_{-0.02}^{+0.02}$		$0$	$-{0.694}_{-0.002}^{+0.026}$		$0$	$-{0.567}_{-0.010}^{+0.025}$		$0$	${0.000}_{-0.013}^{+0.014}$
$\log \xi$	$3.1$	${3.0984}_{-0.0006}^{+0.0008}$		$3.1$	${3.1074}_{-0.0011}^{+0.0006}$		$3.1$	${3.1074}_{-0.0013}^{+0.0013}$		$3.1$	${3.070}_{-0.003}^{+0.003}$
$A_{\mathrm{Fe}}$	$1$	${1.038}_{-0.010}^{+0.009}$		$1$	${1.070}_{-0.009}^{+0.016}$		$1$	${1.085}_{-0.021}^{+0.018}$		$1$	${0.994}_{-0.015}^{+0.032}$
$\Gamma$	$2$	${2.0024}_{-0.0004}^{+0.0003}$		$2$	${1.9965}_{-0.0005}^{+0.0009}$		$2$	${1.9962}_{-0.0005}^{+0.0011}$		$2$	${2.0008}_{-0.0030}^{+0.0005}$
$E_{\mathrm{cut}}$ [keV]	$300$	${300}^{\star }$		$300$	${300}^{\star }$		$300$	${300}^{\star }$		$300$	${300}^{\star }$
$R_f$	–	$>8$		–	$>8$		–	$>9$		–	${0.698}_{-0.042}^{+0.013}$
${\chi }^2/\nu$		$1185.88/1169$			$1092.35/1169$			$1107.52/1169$			$1145.82/1169$
		= 1.01444			= 0.934429			= 0.947411			= 0.980175

Table 2. Table 2: Input parameters and best-fit values for simulations E and F. The reported uncertainties correspond to 90% confidence level for one relevant parameter. ⋆ indicates that the parameter is frozen in the fit. † indicates that the parameter is unconstrained and we just report the best-fit value.

	Simulation E			Simulation F
	Input	Fit		Input	Fit
tbabs
$N_H/{10}^{20}$ cm-2	$6.74$	${6.74}^{\star }$		$6.74$	${6.74}^{\star }$
relxill_nk
$q$	$9$	${5.136}_{-0.035}^{+0.025}$		$9$	${3.847}_{-0.005}^{+0.006}$
$i$ [deg]	$35$	${28.5}_{-2.0}^{+0.5}$		$35$	$\sim 3^{†}$
$a_{\ast }$	$0.95$	${0.9921}_{-0.0012}^{+0.0005}$		$0.95$	${0.9131}_{-0.0007}^{+0.0005}$
${\alpha }_{13}$	$0.3$	$-{0.61}_{-0.06}^{+0.11}$		$-0.3$	$-{1.85}_{-0.03}^{+0.04}$
$\log \xi$	$3.1$	${3.0934}_{-0.0019}^{+0.0016}$		$3.1$	${3.0973}_{-0.0007}^{+0.0007}$
$A_{\mathrm{Fe}}$	$1$	${1.016}_{-0.020}^{+0.024}$		$1$	${1.130}_{-0.010}^{+0.010}$
$\Gamma$	$2$	${2.0072}_{-0.0011}^{+0.0007}$		$2$	${2.0149}_{-0.0003}^{+0.0003}$
$E_{\mathrm{cut}}$ [keV]	$300$	${300}^{\star }$		$300$	${300}^{\star }$
$R_f$	–	$>8$		–	${4.58}_{-0.23}^{+0.04}$
${\chi }^2/\nu$		$1192.34/1169$			$1447.23/1169$
		=1.01997			=1.23801

Equations

$$d{s}^2=g_{tt}d{t}^2+g_{rr}d{r}^2+g_{\theta \theta }d{\theta }^2+2g_{t\varphi }dtd\varphi +g_{\varphi \varphi }d{\varphi }^2,$$

$$T^{\mu \nu }=(\rho +P)u^{\mu }{u}^{\nu }+g^{\mu \nu }P,$$

$$u^{\mu }=(u^t,0,0,u^{\varphi }).$$

$$u_t$$

$$\Omega$$

$$l$$

$${\nabla }_{\nu }{T}^{\mu \nu }=0,$$

$$a^{\mu }=-\frac{g^{\mu \nu }+u^{\mu }{u}^{\nu }}{\rho +P}{\partial }_{\nu }P.$$

$$a_{\mu }={\partial }_{\mu }W,W(P)=-{\int }^P\frac{d{P}^{\prime }}{\rho (P^{\prime })+P^{\prime }}.$$

$$W=W_{\mathrm{in}}+ln\frac{u_t}{u_t^{\mathrm{in}}}+{\int }_{l_{\mathrm{in}}}^l\frac{\Omega d{l}^{\prime }}{1-\Omega l^{\prime }},$$

$$W=\ln (-u_t)+\mathrm{constant}.$$

$$g=\frac{\sqrt{-g_{tt}-2g_{t\varphi }\Omega -g_{\varphi \varphi }{\Omega }^2}}{1+\lambda \Omega },$$

$$N(E_{\mathrm{obs}})=\frac{1}{E_{\mathrm{obs}}}\int g^3{I}_{\mathrm{e}}(E_{\mathrm{e}})\frac{dXdY}{D^2},$$

$$I_{\mathrm{e}}(E_{\mathrm{e}})\propto \frac{\delta (E_{\mathrm{e}}-E_{\mathrm{line}})}{r^q},$$

$$d{s}^2$$

$$\Sigma$$

$$\Delta$$

$$B$$

$$A_1$$

$$A_2$$

$$A_5$$

$$f$$

Full text

††pubyear: 2019††pagerange: Reflection spectra of thick accretion disks–A

1 Introduction

X-ray spectra of accreting black holes often show relativistic reflection features that are thought to be produced by illumination of the inner part of the accretion disk by a hot corona (Fabian et al., 1989; Tanaka et al., 1995; Fabian et al., 2000; Zoghbi et al., 2010; Walton et al., 2012, 2013; Risaliti et al., 2013; Bambi, 2018). Thermal photons from the accretion disk can inverse Compton scatter off free electrons in the corona, which is some hot, usually compact and optically thin, cloud in the strong gravity region near the black hole (e.g. the accretion flow between the disk and the black hole, the atmosphere above the disk, the base of the jet, etc.). A fraction of the Comptonized photons illuminates the accretion disk, generating the reflection component. The latter is characterized by some narrow fluorescent emission lines in the soft X-ray band, in particular the iron K$\alpha$ complex at 6.4 to 6.97 keV depending on the ionization of iron atoms, and by the Compton hump peaked at 20-30 keV. In the presence of high quality data and with the correct astrophysical model, X-ray reflection spectroscopy can be quite a powerful tool to study the strong gravity region of black holes (Brenneman, 2013; Reynolds, 2014) and even to test Einstein’s theory of general relativity (Bambi, 2013; Cao et al., 2018; Abdikamalov et al., 2019b).

Theoretical models to fit these reflection components are typically the combination of a non-relativistic reflection model calculating the reflection spectrum in the rest-frame of the gas, e.g. reflionx (Ross & Fabian, 2005) or xillver (García & Kallman, 2010; García et al., 2013), and a convolution model taking into account the structure of the accretion disk and the relativistic effects of the spacetime, e.g. kerrconv (Brenneman & Reynolds, 2006) or relconv (Dauser et al., 2010, 2013). All the current convolution models employed to fit X-ray data of accreting black holes assume that the disk is infinitesimally thin, perpendicular to the black hole spin, and that the inner edge is at the innermost stable circular orbit (ISCO) or at a larger radius. However, in reality accretion disks around black holes have a finite thickness, which should increase with the mass accretion rate. The inner edge of the disk is thought to be at the ISCO for sources in the high-soft state with an accretion luminosity between 5% to about 30% of their Eddington limit (Steiner et al., 2010; McClintock et al., 2014). For higher mass accretion rate, the gas pressure becomes more and more important, and the disk moves from thin to slim and then thick. The thickness of the disk increases, so the approximation of an infinitesimally thin disk gets worse and worse, and the inner edge of the disk may move inside the ISCO.

While the available relativistic reflection models should only be used to fit sources accreting from thin disks, they are commonly employed to fit any relativistic reflection spectrum (Brenneman, 2013). In particular in the case of supermassive black holes, we know that some sources have accretion luminosities well beyond 30% of their Eddington limit, with some sources even accreting at the Eddington rate and above. This clearly has an impact on the estimation of the model parameters of the system. For example, this may lead to overestimating the black hole spin, because the inner edge of the disk is likely inside the ISCO and the reflection spectrum from the inner part of the disk may look like that of a disk with an inner edge at the ISCO in a spacetime with a smaller ISCO radius (and thus higher black hole spin). The issue is even more subtle if we want to analyze the reflection spectrum of a source to test the Kerr nature of the compact object (Xu et al., 2018; Tripathi et al., 2018, 2019a), because the employment of the wrong model may induce apparent deviations from the Kerr metric (Zhang et al., 2019; Liu et al., 2019).

The aim of the present work is to calculate reflection spectra of thick disks and estimate the impact of the disk structure on the measurement of the parameters of the system when we fit the data with a thin disk model. In particular, we are interested in the impact on the estimation of the deformation parameters of the spacetime when we want to test the Kerr metric. Surprisingly, there is not much on this subject in the literature, despite it being quite clear that this effect exists and can lead to some erroneous estimates of the properties of accreting black holes. Wu & Wang (2007) present iron line profiles for thick and sub-Keplerian accretion disks with inner edges at the ISCO or at larger radii. Taylor & Reynolds (2018) study the reflection spectrum from thin disks with inner edges at the ISCO radius for a disk model with finite thickness. In our study, we describe the accretion disk with the Polish donut model (Kozlowski et al., 1978; Abramowicz et al., 1978), which is expected to be suitable for objects with very high accretion rates, say around the Eddington limit. Here the thickness of the disk is $H/R\sim 1$, where $H$ is the semi-height of the disk at the radius $R$. The inner edge of the disk is always inside the ISCO and can approach the marginally bound orbit.

We find that single line shapes of thick disks are definitively different from those generated by thin disks. However, when we consider the full reflection spectrum and we do not know the input values of the model parameters, it is very challenging to recover the correct values of the model parameters from the fit. In other words, some parameters show strong degeneracy, and the reflection spectrum of a thick disk can look like that of a thin disk with different values of some model parameters. Since the inner edge of the disk is inside the ISCO radius in the case of Polish donut disks, the spin parameter is usually overestimated. In the case of tests of the Kerr metric, a Kerr black hole may be interpreted as a non-Kerr compact object. However, as the spin parameter increases, real deviations from the Kerr background have quite a strong impact on the reflection spectrum and, when we fit the spectrum of a thick disk of a Kerr black hole with a thin disk model, we may not measure deviations from the Kerr geometry even if we still overestimate the black hole spin.

The content of the paper is as follows. In Section 2, we briefly review the Polish donut model, which is the model adopted in this work to describe thick accretion disks. In Section 3, we present our convolution model and we show iron line shapes expected from thick disks described by the Polish donut model for different values of the black hole spin and of the inclination angle of the disk. In Section 4, we consider the whole reflection spectrum and we simulate some observations with NICER (Gendreau et al., 2016) to determine the impact of the disk structure on the estimation of the black hole spin and the deformation parameter of the spacetime. We discuss our results in Section 5. In Appendix A, we review the Johannsen metric (Johannsen, 2013), which is the spacetime employed in the reflection model used in Section 4 to analyze the simulations and estimate the impact of the disk structure on the measurement of the model parameters. Throughout the manuscript, we adopt units in which $c=G_{\rm N}=1$ and a metric with signature $(-+++)$.

2 Accretion disk model

Analytic models for relativistic thick accretion disks with negligible self-gravitation were originally developed in Fishbone & Moncrief (1976), Kozlowski et al. (1978), and Abramowicz et al. (1978). These authors considered stationary, axially symmetric, perfect fluid disks in generic stationary, axisymmetric, and asymptotically flat spacetimes. However, Fishbone & Moncrief (1976) assumed isentropic disks, while Kozlowski et al. (1978) and Abramowicz et al. (1978) worked with barotropic disks. The latter is today known as the “Polish donut” model and describes a low viscosity, geometrically and optically thick, radiation pressure supported disk undergoing super-Eddington black hole accretion. In this section, we just provide a brief overview. More details can be found, for instance, in Li & Bambi (2013) and Bambi (2017b).

We assume a stationary, axisymmetric, and asymptotically flat spacetime. We write the line element as111We note that the most general stationary, axisymmetric, and asymptotically flat spacetime also has a non-vanishing $g_{tr}$ [see, for instance, Section 7.1 in Wald (1984)], but here we ignore such a possibility. In general relativity in vacuum, we can always find a coordinate system with the line element as in Eq. (1) where $g_{tr}=0$.

[TABLE]

where the metric elements are all independent of the $t$ and $\phi$ coordinates. The matter in the disk is treated as a perfect fluid with purely azimuthal flow. The stress-energy tensor of the system can be written as

[TABLE]

where $\rho$ and $P$ are, respectively, the energy density and the pressure of the perfect fluid, which are related through a barotropic equation of state $\rho=P(\rho)$. The 4-velocity vector of the fluid element is

[TABLE]

The specific energy of the fluid element is $-u_{t}$. The angular velocity and the angular momentum per unit energy are $\Omega=u^{\phi}/u^{t}$ and $l=-u_{\phi}/u_{t}$, respectively. $u_{t}$, $\Omega$, and $l$ can be written in terms of metric elements as

[TABLE]

The motion of a fluid element in the disk is governed by Euler’s equation

[TABLE]

which leads to the following equation for the 4-acceleration of the fluid element $a^{\mu}=u^{\nu}\nabla_{\nu}u^{\mu}$

[TABLE]

Exploiting the assumptions that the spacetime is stationary and axisymmetric and that the equation of state is barotropic, the 4-acceleration of the fluid element can be written as a gradient of the scalar potential $W(P)$

[TABLE]

From Eq. (7) we can infer that there exists an invariant function $\Omega=\Omega(l)$ (Abramowicz et al., 1978). Integrating Eq. (7), we find

[TABLE]

where $W_{\rm in}$, $-u_{t}^{\rm in}$ and $l_{\rm in}$ are, respectively, the potential, the energy per unit mass, and the angular momentum per unit energy at the inner edge of the disk and can be replaced by the same quantities at any other point of the fluid boundary. In the Newtonian limit, $W$ reduces to the sum of the gravitational and the centrifugal potential and vanishes at infinity. For the special case $l={\rm constant}$, Eq. (10) still holds, and the integral vanishes.

Once the background metric is known, the function $\Omega=\Omega(l)$ determines the dissipative process and characterizes the fluid rotation in the disk as well. In the zero viscosity case, this function cannot be determined from any equation and it should come from the assumption on the model. If we fix $\Omega=\Omega(l)$, we can determine the surfaces of the constant pressure of the fluid, which are usually called the equipotential surfaces and have $W={\rm constant}<0$. Such surfaces describe possible boundaries of the fluid. Accretion is possible when the fluid surface has one (or more) cusp(s): the fluid fills out the Roche lobe and is then transferred to the compact object (Abramowicz et al., 1978; Li & Bambi, 2013; Bambi, 2017b). The mechanism does not need any fluid viscosity, which, at least in principle, could thus be very low.

For simplicity, in this work we consider the special case $l={\rm constant}$, which is marginally stable with respect to axisymmetric perturbations. For $l={\rm constant}$, Eq. (10) reduces to

[TABLE]

In the Kerr spacetime, there are five qualitatively different configurations, depending upon the value of $l$. The most interesting and relevant to the present work is the one in which $l_{\rm ms}<l<l_{\rm mb}$, where $l_{\rm mb}$ and $l_{\rm ms}$ are, respectively, the angular momentum per unit energy at the marginally bound and the marginally stable (or ISCO) orbits. For these values of $l$, there is only one surface with cusp $(W=W_{\rm cusp})$ and the location of the cusp lies in between the marginally stable and marginally bound orbits. The configurations with surface $W<W_{\rm cusp}$ describe non-accreting disks, while those with surface $W>W_{\rm cusp}$ do not represent any disk configuration (Abramowicz et al., 1978; Li & Bambi, 2013; Bambi, 2017b).

3 Iron line shapes

Let us consider a black hole accreting from a geometrically and optically thick disk, which we will describe with the Polish donut model. The illumination of the disk by the corona produces a reflection component. In order to understand the impact of the disk structure on the reflection spectrum, in this section we will simplify the whole reflection spectrum with a single iron line at 6.4 keV. Full reflection spectra will be considered in the next section.

Calculations of single iron lines of thin accretion disks have already been described in the literature; see, e.g., Bambi (2013) or Bambi (2017b). We consider a distant observer with viewing angle $i$ and we calculate backward in time photon trajectories from the image plane of the distant observer to the emission points on the accretion disk. When a photon hits the disk, we calculate the redshift factor

[TABLE]

where $\Omega$ is the angular velocity of the fluid element, $\lambda=k_{\phi}/k_{t}$, and $k_{t}$ and $k_{\phi}$ are, respectively, the $t$ and $\phi$ components of the photon 4-momentum. In a stationary and axisymmetric spacetime, $k_{t}$ and $k_{\phi}$ are constants of motion, and therefore $\lambda$ is a constant of motion too and can be evaluated from the photon initial conditions. Finally, the iron line shape of the accretion disk is obtained by integrating over the disk image

[TABLE]

where $N(E_{\rm obs})$ is the photon count with energy $E_{\rm obs}$ (as measured by the distant observer), $I_{\rm e}$ and $E_{\rm e}$ are, respectively, the specific intensity of radiation and the photon energy at the emission point on the disk (as measured in the rest frame of the accreting gas), $X$ and $Y$ are the Cartesian coordinates of the image plane of the distant observer, and $D$ is the distance of the source from the distant observer. In the case of monochromatic emission (single line) with power-law profile, the specific intensity of radiation $I_{\rm e}$ can be written as

[TABLE]

where $E_{\rm line}=6.4$ keV for a neutral iron K$\alpha$ line and $q$ is the emissivity index. More details can be found, for instance, in Bambi (2013) or Bambi (2017b).

In the case of a thick disk described by the Polish donut model, the calculations are quite similar. We still start from the image plane of the distant observer and we calculate backward in time the null geodesics of a grid of photons to obtain the image map of the accretion disk, with photon redshift and photon intensity at every pixel of the image. For a specific metric, we fix the value of $l$ and we infer the boundary of the accretion disk. Photon trajectories are still calculated until the photon hits the disk, but now the emission point is outside the equatorial plane and we need a more sophisticated routine to determine its coordinates. Once the coordinates of the emission point are determined, we still employ Eq. (12) to infer the redshift factor. Now $\Omega$ is calculated from Eq. (5), while in the thin disk case $\Omega$ was the Keplerian angular velocity. Note also that some parts of the disk may not be observed as they are obscured by the disk itself; the phenomenon does not occur for face-on disks, becomes more and more important as the viewing angle increases, and the inner part of the disk is completely obscured for high viewing angles.

The ray-tracing code to calculate the photon trajectories from the image plane of the distant observer to the emission point on the accretion disk is that described in Abdikamalov et al. (2019a). However, unlike the code in Abdikamalov et al. (2019a) that is used to calculate the transfer function of the spacetime, we directly evaluate the photon count per energy bin using Eq. (14), as done in Bambi (2012). In such a case, it is easier to reach a better accuracy since we skip the numerical uncertainties related to the calculation of the transfer function, its tabulation, and its integration. We have compared our calculations for thin disks with relline (Dauser et al., 2010, 2013), finding that the discrepancy in the photon flux is below 0.1%, see Fig. 2. We cannot do the same for our calculations for thick disks because there are no existing models for a similar check.

Iron line shapes calculated by our ray-tracing code for $q=3$ and 9 are shown, respectively, in Fig. 3 and Fig. 4 for a dimensionless spin parameter $a_{*}=0$, 0.6, 0.9, 0.95, and 0.998 (from top to bottom) and for a viewing angle $i=10^{\circ}$, $35^{\circ}$, and $60^{\circ}$ (from left to right). In every panel, the black dotted-dashed curve indicates the iron line of a thin disk with inner edge at the ISCO and outer edge $R_{\rm out}=400$ $M$. For the Polish donut disks, we do not fix $l$ but the outer edge of the disk $R_{\rm out}$ in order to compare similar size disks in different spacetimes. The outer edge of the disk is at $R_{\rm out}=12$ $M$ (blue solid curves), $20$ $M$ (green dotted curves), and $40$ $M$ (orange dashed curves). Such values of the outer edges are unnaturally low for real accretion disks, but, within the Polish donut model, larger disks do not permit the observation of their inner part except for a viewing angle $i$ close to $0^{\circ}$ (face-on disks), see Fig. 2. This fact challenges the comparison between our thin and thick disks because of the finite size effects of the thick disks, but larger thick disks would make the comparison completely impossible because their spectrum does not include the strongly redshifted photons from the region of the disk closer to the black hole. The choice of the high emissivity index $q=9$ is to minimize such a problem related to the Polish donut model, as with such a high value of $q$ most of the radiation is emitted from the inner edge. The inner edge of the disk is always inside the ISCO. In Fig. 3 and Fig. 4 we see that the iron lines of the polish donut model have some bumpiness, particularly at high inclination. This also shows up in the thin disk lines at higher spin in Fig. 4. This is just the result of some numerical uncertainty in the calculations, but it is completely negligible for what follows.

The iron line shapes of thick accretion disks are remarkably different from those of thin disks in the same spacetimes and with the same inclination angles. For $q=3$, the intensity of the radiation decays slower as the radial coordinate increases, so the effects of the outer edge of the disk in the Polish donut model are more pronounced. For $q=9$, most of the radiation comes from the very inner part of the accretion disk, so the differences in the line shapes is mainly determined by disk structure near the inner edge. For thicker disks (larger $R_{\rm out}$) and/or higher viewing angles, the inner part of the disk may be partially or completely obscured, with the result that the low energy tail of the line (due to strongly redshifted photons emitted near the inner edge of the disk) is missing. This can be quite clearly seen from the orange curves for $i=60^{\circ}$, both in Fig. 3 and Fig. 4. If we increase the value of the spin parameter $a_{*}$, the inner edge of the disk moves to smaller radii. In general, this should produce highly redshifted photons. However, in the case of thick disks such photons may not reach the distant observer, except in the case of almost face-on disks, so a higher spin parameter does not necessarily imply a more redshifted iron line.

4 Reflection spectra, simulations, and impact on parameter estimates

All the available relativistic reflection models assume infinitesimally thin accretion disks with the inner edge at the ISCO or at a larger radius. In principle, such models should only be used to fit the reflection spectra of sources with thin disks. In reality, these models are commonly used to fit relativistic reflection features of any source, regardless of its mass accretion rate. This, of course, leads to systematic uncertainties in the final measurements of the model parameters. For instance, Tab. 1 in Brenneman (2013) shows a summary of black hole spin measurements derived from the analysis of reflection features and an estimate of the accretion luminosity of the corresponding source222We note that Tab. 1 in Brenneman (2013) shows the total luminosity of every source. A better proxy for the accretion luminosity and the disk thickness would be the disk luminosity.. If the latter are correct, some accretion disks must be definitively thick.

In our case, we are particularly interested in the impact of the disk structure on tests of the Kerr metric (Bambi, 2017a). In a certain sense, it is indeed surprising that it is possible to obtain very strong constraints on possible deviations from the Kerr metric using thin disk relativistic reflection models from supermassive black holes that should not accrete from thin disks, see Tripathi et al. (2019a). Naively, we could expect that the systematic uncertainties of the model can lead to measurements of deviations from the Kerr spacetime. On the contrary, in Tripathi et al. (2019a), we intentionally selected very-fast rotating black holes (with no preference based on accretion luminosity and some sources are indeed thought to exceed their Eddington limit) and we found very stringent constraints on possible deviations from the Kerr background. Here, we want to understand how this is possible by simulating some reflection spectra of thick disks and by fitting the faked data with our thin disk relativistic reflection model relxill_nk to test the Kerr hypothesis (Bambi et al., 2017; Abdikamalov et al., 2019a)333relxill_nk is an extension of the relxill package (Dauser et al., 2013; García et al., 2013) to non-Kerr spacetimes and the two last letters NK stands for Non-Kerr. See Bambi et al. (2017); Abdikamalov et al. (2019a) for more details..

To simulate reflection spectra, we use the ray-tracing code described in the previous section to calculate the photon flux via Eq. (13). Now the specific intensity of the radiation at the emission point, $I_{\rm e}$, is not given by Eq. (14) but is calculated with xillver (García & Kallman, 2010; García et al., 2013)444xillver is a non-relativistic reflection model: it calculates the reflection spectrum generated by the illumination of a cold material (not necessarily a disk) by an incident radiation. Disk geometry and spacetime metric are included after convolving the output of xillver with a convolution model. In the Kerr metric for an infinitesimally thin disk, this is done with relconv and the result is relxill: in XSPEC language, relxill = relconv $\times$ xillver. In this work, we simply replace relconv with a convolution model for thick disks.. Once we have the theoretical reflection spectrum of a given configuration (specified by the parameters of the spacetime, of the disk, of xillver, and by the viewing angle), we use XSPEC and the fakeit command to generate the simulated observations. For the X-ray mission, we choose NICER (Gendreau et al., 2016), and therefore we use its background, ancillary, and response matrix files. NICER has a good energy resolution near the iron line, which make it suitable for X-ray reflection spectroscopy. It was launched in 2017 and there is already some important work in literature on spin measurements of black hole binaries using the iron line method with this mission.

In XSPEC language, the simulations are done with the model

tbabs $\times$ ( powerlw + reflection ) ,

where tbabs takes the Galactic absorption into account (Wilms et al., 2000), powerlw describes the continuum from the corona, and reflection is the reflection component obtained by applying our convolution model to the output of xillver.

First, we consider four simulations in the Kerr spacetime. The black hole spin parameter, $a_{*}$, is set to 0.8 (simulation A), 0.85 (simulation B), 0.9 (simulation C), and 0.95 (simulation D). For all the simulations, the inclination angle is fixed to $35^{\circ}$ and in xillver we adopt the ionization $\log\xi=3.1$ ($\xi$ in units erg cm s*-1*), the iron abundance $A_{\rm Fe}=1$ (Solar abundance), and the photon index $\Gamma=2$. We set the outer edge of the accretion disk at $R_{\rm out}=40$ $M$ for all simulations. The inner edges of the disk turns out to be $R_{\rm in}=2.129$ $M$ (simulation A), 1.951 $M$ (simulation B), 1.750 $M$ (simulation C), and 1.507 $M$ (simulation D). The emissivity profile of the reflection spectrum is described by a power-law with emissivity index $q=9$ in order to limit the effects of the outer edge of the disk. We set the photon flux around $1.4\cdot 10^{-10}$ erg cm*-2* s*-1* and the exposure time around 420 ks in order to have about 50 million photons in the energy range 0.2-12 keV. Such a photon count is higher than the typical photon count for an observation of a bright active galactic nucleus, so we can potentially better constrain the model parameters than what it is possible today with current X-ray observations.

The simulated data are analyzed with the XSPEC model

tbabs $\times$ relxill_nk ,

where relxill_nk is our relativistic reflection model (Bambi et al., 2017; Abdikamalov et al., 2019a). We do not assume the Kerr spacetime and we employ the Johannsen metric with a possible non-vanishing deformation parameter $\alpha_{13}$ (Johannsen, 2013). The Johannsen metric is briefly reviewed in Appendix A. It is not a solution of a well-defined theory of gravity but a parametric black hole metric in which an infinite number of deformation parameters are introduced to quantify possible deviations from the Kerr background. When all deformation parameters vanish, we exactly recover the Kerr solution. With such a metric, it is like performing a null experiment: we want to test the Kerr nature of a black hole and we try to measure the deformation parameters of the Johannsen metric to check whether they vanish.

Tab. 1 shows the results of our analysis of simulations A-D. The quality of the best-fit can be seen from Fig. 6, which shows the ratios between the data and the best-fit models. The results of the analysis of these simulations are discussed in the next section.

We consider two more simulations imposing a non-Kerr background. The black hole spin parameter is set to 0.95 and the deformation parameter is $\alpha_{13}=0.3$ (simulation E) and $-0.3$ (simulation F). The values of the other model parameters are the same as simulations A-D. Even for simulations E and F we use NICER and we require about 50 million photons in the energy range 0.2-12 keV. We analyze the data with relxill_nk and we report best-fit values and ratio plots in Tab. 2 and in the two bottom panels of Fig. 6, respectively.

4.1 Slab reflector vs Toroidal reprocessor

For the calculation of reflection spectra, we have assumed the Polish donut is optically thick through the use of a cold slab reflector. However, there can be scenarios in which the disk is not optically thick. For this reason, here we perform a brief comparison of the iron line generated by the optically thick model we have used throughout this work to that of a model that includes the effect of scattering that would be present in a disk that is not totally opaque. For the scattering model we use the Compton-thick toroidal reprocessor created by Murphy & Yaqoob (2009) called MYTorus. The geometry used is that of a torus with a circular cross-section and is placed such that it subtends a solid angle of $2\pi$ at the central source, which emits isotropically. The model includes the effect of both the reflected and scattered emission from the cold (neutral) reprocessor (for more details see Murphy & Yaqoob (2009)). Since the X-ray source used in MYTorus is a centrally-located point source, while in the model used throughout this work nothing is assumed about the location or geometry of the source the comparison here is not exact but rather an approximate look at the effect of including scattering.

For the comparison we will only study the effect on the iron line as the difference is much more noticeable there than in the full reflection spectrum. We compare the following four models in Fig. 5

MYTorus iron line 2. 2.

Polish donut iron line 3. 3.

relconv $\times$ MYTorus iron line 4. 4.

Polish donut $\times$ MYTorus iron line

where "iron line" means we only take the iron line part of the model and $\times$ refers to a convolution. The blue solid lines are the Fe K$\alpha$ and K$\beta$ emission lines produced by the model MYTorus. These lines are intrinsically narrow because of the fact that MYTorus model does not include relativistic effects. These lines contain both the zeroth-order photons as well as the scattered photons. The input parameters used to generate these lines are: column density $N_{H}=2\times 10^{24}$ cm*-2*, photon index $\Gamma=2$, redshift $z=0$. We used the version of the fits table to generate the iron lines for which offset energy of the line is 0 and the termination energy of the incident power-law continuum is 300 keV. In order to introduce the relativistic effects in the MYTorus iron lines we convolved it with the convolution model relconv. The parameters used for relconv are: spin $a*=0.95$, emissivity index $q=9$ and outer edge of the disk $R_{\rm out}=40$ $M$. In this case the relativistic effects are more important than the geometry of the disk. The resultant iron line becomes broadened (green dashed curve, Fig. 5). The effects of convolving the iron lines produced by MYTorus with the Polish donut model are shown by the orange dotted-dashed curves. In this case the geometry of the disk becomes important and the iron lines show significant differences as compared to the relconv$\times$ MYTorus iron line (green dashed curves). However, the differences are negligible when the Polish donut $\times$ MYTorus iron line is compared with the iron lines produced by the Polish donut model. Thus, it seems that including scattering has some effect on the iron line, and in turn the full reflection spectrum, however this effect is likely not large enough to have a significant impact on recovered parameters.

5 Discussion and conclusions

In the previous section, we have simulated observations of reflection spectra from thick accretion disks in Kerr (simulations A-D) and non-Kerr (simulations E and F) spacetimes and we have analyzed the data with the thin disk relativistic reflection model relxill_nk to figure out potential systematic uncertainties when we fit the reflection spectrum of a source accreting at Eddington or super-Eddington rate with a thin disk model. Note that this is the case shown in many studies in the literature of supermassive black holes. The Polish donut disk employed in our simulations is quite a simple model far from a real accretion disk, but it is used here to catch the basic differences between thin and thick disks and the analytic model helps to simplify our study.

Let us first discuss the results of simulations A-D. From Tab. 1 and the ratio plots in Fig. 6, we can say:

The quality of the fits is good: the reduced $\chi^{2}$ is close to 1 and we do not see unresolved features in the ratio plots. In other words, reflection spectra from thick disks can look like those from thin disks. This may sound a little bit surprising at first, because from Figs. 3 and 4 we see that the iron line shapes of thick disks are very different from the same lines generated by thin disks. However, in Figs. 3 and 4 all the other model parameters are fixed and we only change the disk structure. In reality, we have several free parameters and, eventually, reflection spectra from thick disks can look like those from thin disks with different model parameters. 2. 2.

The fit can recover correct values of some model parameters and wrong values of other model parameters. For example, the fit can measure more or less the right value of $\log\xi$, $A_{\rm Fe}$, and $\Gamma$ even if we employ a thin disk model to fit the spectrum of a thick disk. On the contrary, the estimates of the emissivity index $q$, the spin parameter $a_{*}$, and the deformation parameter $\alpha_{13}$ are strongly affected by the systematic uncertainties related to the use of a wrong model. 3. 3.

The spin parameter $a_{*}$ is always and significantly overestimated. This result could be expected, because the inner edge of the disk is inside the ISCO in all the simulations and therefore some photons are strongly redshifted, as it would be the case of emission from the ISCO in a spacetime with smaller ISCO radius. 4. 4.

While the simulations are done in a Kerr spacetime, the fit “typically” measures a non-vanishing value of the deformation parameter $\alpha_{13}$. Generally speaking, this is also a result that could be expected: we are fitting the data with a wrong model and therefore we find a wrong estimate of the deformation parameter. Interestingly, the best-fit of $\alpha_{13}$ is far from the Kerr solution for simulation A and moves to the Kerr solution $\alpha_{13}=0$ as the spin parameter increases, to recover the Kerr metric in simulation D (but with too large a value of the spin parameter with respect to its input value in the simulation). We are tempted to argue that, when the inner edge of the disk is very close to the black hole, relativistic effects in the reflection spectrum are very important so that even small deviations from Kerr could produce very different reflection spectra. Note also that the results of simulation D seem to nicely explain the very strong constraints on the Kerr hypothesis found in Tripathi et al. (2019a), where we analyzed some bare active galactic nuclei: some sources look like near-extremal Kerr black holes, with the best-fit of the spin parameter $a_{*}$ stuck at the maximum value of relxill_nk, namely 0.998, and very good constraints on the deformation parameter $\alpha_{13}$. This is exactly what we find in simulation D. If such an interpretation is correct, we should be very skeptical of the extremely high values of the spin parameters found in Tripathi et al. (2019a) for those sources.

Simulations E and F in non-Kerr spacetime mainly confirm the conclusions derived from the results of simulations A-D. The fit can recover the correct value of some parameters ($\log\xi$, $A_{\rm Fe}$, and $\Gamma$) and infers wrong values for other parameters ($q$, $a_{*}$, and $\alpha_{13}$). For simulations F, in which the input value is $\alpha_{13}=-0.3$, the quality of the fit is not good: the reduced $\chi^{2}$ is around 1.2 and its ratio plot shows some unresolved features. In all simulations, in Kerr and non-Kerr spacetimes but with the exception of simulation D, the best-fit value of $\alpha_{13}$ is negative. This, presumably, is just a property of the deformation parameter $\alpha_{13}$, which can somehow mimic the effect of a Polish donut disk on the reflection spectrum when it becomes negative. For simulation D, the inner edge of the disk is probably so close to the black hole event horizon that relativistic effects are more important than the disk structure. This would explain the trend shown in Tab. 1, where the estimate of $\alpha_{13}$ moves to the Kerr case $\alpha_{13}=0$ as the inner edge of the disk gets closer to the black hole event horizon.

In conclusion, the results of our work suggest that black hole spin measurements and tests of the Kerr hypothesis from the study of relativistically broadened reflection features can be significantly affected by systematic uncertainties when we employ a thin disk model to analyze sources accreting at the Eddington or super-Eddington rate, as is often the case for many supermassive black holes.

Acknowledgements

This work was supported by the Innovation Program of the Shanghai Municipal Education Commission, Grant No. 2019-01-07-00-07-E00035, and Fudan University, Grant No. IDH1512060. S.N. acknowledges support from the Alexander von Humboldt Foundation and the Excellence Initiative at Eberhard-Karls Universität Tübingen.

Appendix A Johannsen metric

The Johannsen background is a parametric black hole spacetime designed to quantify possible deviations from the Kerr solution (Johannsen, 2013). In Boyer-Lindquist-like coordinates, its line element reads

[TABLE]

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

[TABLE]

The functions $A_{1}$, $A_{2}$, $A_{5}$, and $f$ are defined as

[TABLE]

and $\{\alpha_{1n}\}$, $\{\alpha_{2n}\}$, $\{\alpha_{5n}\}$, and $\{\epsilon_{n}\}$ are four sets of deformation parameters without constraints from the Newtonian limit and experiments in the Solar System. In the present manuscript, as a preliminary study on the impact of the disk structure on tests of the Kerr hypothesis, we have only considered the possibility of a non-vanishing $\alpha_{13}$ and set to zero all other deformation parameters.