A CODE TO COMPUTE THE EMISSION OF THIN ACCRETION DISKS IN NON-KERR SPACETIMES AND TEST THE NATURE OF BLACK HOLE CANDIDATES

The calculation of the thermal spectrum of a thin accretion disk has been extensively discussed in the literature, see, e.g., Li et al. (2005) and references therein. The spectrum can be conveniently written in terms of the photon flux number density as measured by a distant observer, $N_{E_{\rm obs}}$. A quite common approximation in the literature is to neglect the effect of light bending, which is good for small inclination angles of the disk (face-on disks). In this case, $N_{E_{\rm obs}}$ is

$$\begin{eqnarray} N_{E_{\rm obs}} &=& \frac{1}{E_{\rm obs}} \int I_{\rm obs}(\nu) d \Omega _{\rm obs} \nonumber \\ &=& \frac{1}{E_{\rm obs}} \int g^3 I_{\rm e}(\nu _{\rm e}) d \Omega _{\rm obs} \nonumber \\ &=& A_1 \left(\frac{E_{\rm obs}}{\rm keV}\right)^2 \cos i \int _{\rm disk} \frac{1}{M^2}\frac{\Upsilon \sqrt{-G} dr d\phi }{\exp \left[\frac{A_2}{g F^{1/4}} \left(\frac{E_{\rm obs}}{\rm keV}\right)\right] - 1},\,\, \quad \end{eqnarray} \tag{ 4 }$$

where I_obs, E_obs, and ν are, respectively, the specific intensity of the radiation, the photon energy, and the photon frequency measured by the distant observer. $d\Omega _{\rm obs} \approx \cos i \sqrt{-G} dr d\phi / D^2$ is the element of the solid angle subtended by the image of the disk on the observer's sky, i is the viewing angle of the observer, and D is the distance of the source. g is the redshift factor

$$\begin{eqnarray} &&g = \frac{E_{\rm obs}}{E_{\rm e}} = \frac{\nu }{\nu _{\rm e}} = \frac{k_\alpha u^{\alpha }_{\rm obs}}{k_\beta u^{\beta }_{\rm e}}, \end{eqnarray} \tag{ 5 }$$

where E_e = hν_e, k^α is the 4-momentum of the photon, u^α_obs = (− 1, 0, 0, 0) is the 4-velocity of the distant observer, and u^α_e = (u^t_e, 0, 0, Ωu^t_e) is the 4-velocity of the emitter. I_e(ν_e)/ν³_e = I_obs(ν_obs)/ν³ follows from Liouville's theorem. A₁ and A₂ are given by (for the sake of clarity, here I show explicitly G_N and c)

$$\begin{eqnarray} A_1 &=& \frac{2 \left({\rm keV}\right)^2}{f_{\rm col}^4} \left(\frac{G_{\rm N} M}{c^3 h D}\right)^2 \nonumber \\ &=& \frac{0.07205}{f_{\rm col}^4} \left(\frac{M}{M_\odot }\right)^2 \left(\frac{\rm kpc}{D}\right)^2 \, {\rm \gamma \, keV^{-1} \, cm^{-2} \, s^{-1}}, \nonumber \\ A_2 &=& \left(\frac{\rm keV}{k_{\rm B} f_{\rm col}}\right) \left(\frac{G_{\rm N} M}{c^3}\right)^{1/2} \left(\frac{4 \pi \sigma }{\dot{M}}\right)^{1/4} \nonumber \\ &=& \frac{0.1331}{f_{\rm col}} \left(\frac{\rm 10^{18} \, g \, s^{-1}}{\dot{M}}\right)^{1/4} \left(\frac{M}{M_\odot }\right)^{1/2}. \end{eqnarray} \tag{ 6 }$$

Using the normalization condition g_μνu^μ_eu^ν_e = −1, one finds

$$\begin{eqnarray} &&u^t_{\rm e} = - \frac{1}{\sqrt{-g_{tt} - 2 g_{t\phi } \Omega - g_{\phi \phi } \Omega ^2}}, \end{eqnarray} \tag{ 7 }$$

and therefore

$$\begin{eqnarray} &&g = \frac{\sqrt{-g_{tt} - 2 g_{t\phi } \Omega - g_{\phi \phi } \Omega ^2}}{1 + \lambda \Omega }, \end{eqnarray} \tag{ 8 }$$

where λ = k_ϕ/k_t is a constant of the motion along the photon path. Neglecting the effect of light bending, λ = rsin ϕsin i, where r and ϕ are the coordinates on the equatorial plane of the emitter. Doppler boosting, gravitational redshift, and frame dragging are entirely encoded in the redshift factor g.

The effect of light bending can be included by using a ray-tracing approach. The initial conditions (t₀, r₀, θ₀, ϕ₀) for the photon with Cartesian coordinates (X, Y) on the image plane of the distant observer are given by (Johannsen & Psaltis 2010)

$$\begin{eqnarray} t_0 &=& 0, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} r_0 &=& \sqrt{X^2 + Y^2 + D^2}, \, \, \, \qquad \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \theta _0 &=& \arccos \frac{Y \sin i + D \cos i}{\sqrt{X^2 + Y^2 + D^2}}, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \phi _0 &=& \arctan \frac{X}{D \sin i - Y \cos i}. \end{eqnarray} \tag{ 12 }$$

As the initial 3-momentum $\bf {k}_0$ must be perpendicular to the plane of the image of the observer, the initial conditions for the 4-momentum of the photon are (Johannsen & Psaltis 2010)

$$\begin{eqnarray} k^r_0 &=& - \frac{D}{\sqrt{X^2 + Y^2 + D^2}} |\bf {k}_0|, \qquad \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \, \, \, k^\theta _0 &=& \frac{\cos i - D \frac{Y \sin i + D \cos i}{X^2 + Y^2 + D^2}}{\sqrt{X^2 + (D \sin i - Y \cos i)^2}} |\bf {k}_0|, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} k^\phi _0 &=& \frac{X \sin i}{X^2 + (D \sin i - Y \cos i)^2} |\bf {k}_0|, \, \, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \quad \, \, k^t_0 &=& \sqrt{\big(k^r_0 \big)^2 + r^2_0 \big(k^\theta _0 \big)^2 + r_0^2 \sin ^2\theta _0 \big(k^\phi _0 \big)^2}. \end{eqnarray} \tag{ 16 }$$

In the numerical calculation, the observer is located at D = 10⁶ M, which is far enough to assume that the background geometry is flat and therefore k^t₀ can be inferred from the condition g_μνk^μk^ν = 0 with the metric tensor of a flat spacetime. The photon trajectory is numerically integrated backward in time to the point of the photon emission on the accretion disk: in this way, we get the radial coordinate r_e at which the photon was emitted and the angle ξ between the wavevector of the photon and the normal of the disk surface (necessary to compute ϒ). Now g = g(X, Y), as in Equation (8) everything depends on r_e, except λ that can be evaluated from the photon initial conditions (λ = k_ϕ/k_t = r₀sin θ₀ k^ϕ₀/k^t₀). The observer's sky is divided into a number of small elements and the ray-tracing procedure provides the observed flux density from each element; summing up all the elements, we get the total observed flux density of the disk. In the case of Kerr background, one can actually exploit the special properties of the Kerr solution and solve a simplified set of differential equations. That is not possible in a generic non-Kerr background and so the code solves the second-order photon geodesic equations of the spacetime, by using the fourth-order Runge–Kutta–Nyström method (Lund et al. 2009). The photon flux number density is given by

$$\begin{eqnarray} N_{E_{\rm obs}} &=& \frac{1}{E_{\rm obs}} \int I_{\rm obs}(\nu) d \Omega _{\rm obs} = \frac{1}{E_{\rm obs}} \int g^3 I_{\rm e}(\nu _{\rm e}) d \Omega _{\rm obs} \nonumber \\ &=& A_1 \left(\frac{E_{\rm obs}}{\rm keV}\right)^2 \int \frac{1}{M^2} \frac{\Upsilon dXdY}{\exp \left[\frac{A_2}{g F^{1/4}} \left(\frac{E_{\rm obs}}{\rm keV}\right)\right] - 1}, \end{eqnarray} \tag{ 17 }$$

where X and Y are the coordinates of the position of the photon on the sky, as seen by the distant observer; that is, dΩ_obs = dXdY/D². The basic features of the code are outlined in Table 1.

Table 1. Basic Features of the Code to Compute the Observed Spectrum of the Disk

Geometry of the background	Johannsen–Psaltis spacetime with three parameters (mass M, spin parameter a, deformation parameter 3). When 3 = 0, the geometry of the spacetime reduces exactly to the Kerr solution. No restrictions on the values of a/M and 3
Relativistic effects	All relativistic effects are included. Ray-tracing technique used
Self-irradiation	Not included
Non-zero torque at rin	Not included
Indirect images	Not included
Color factor fcol	Constant. Set by the user
Radiation emission ϒ	Isotropic or limb-darkened

Note. Indirect images are the images formed by null geodesics penetrating the equatorial plane inside r_in.

Download table as:  ASCIITypeset image

The results of the code in the Kerr background are summarized in Figures 1 and 2, while Table 2 shows some fundamental properties of these spacetimes. The two panels of Figure 1 should be compared with the ones in Figure 5 of Li et al. (2005). The agreement between the two sets of spectra is very good, except for the case a/M = 0.999 and inclination angle i = 30°, for which the difference is however small. The discrepancy is due to the fact that the present code does not include the effect of self-irradiation of the disk. As discussed in Li et al. (2005), the effect of self-irradiation of the disk, as well as a possible non-zero torque at the inner edge of the disk (also not included here, see Table 1), can be ignored in the analysis of observational data, because they can be absorbed by adjusting the mass accretion rate and the spectral color factor of a zero torque model without returning radiation. This is indeed our case, as shown in the left panel of Figure 2, where we can see the spectra of a Kerr BH with a/M = 0.999 and i = 30°, respectively for $\dot{M} = 1.0 \times 10^{19}$ and 1.2 × 10¹⁹ g s⁻¹. The spectrum with higher mass accretion rate and without self-irradiation of the disk is in agreement with the one calculated in Li et al. (2005) for $\dot{M} = 1.0 \times 10^{19}$ g s⁻¹ and with self-irradiation.² The effect of light bending can be seen in the right panel in Figure 2, which shows the observed spectra as computed respectively from Equations (4) and (17). The difference in the observed spectra is significant for large inclination angles (edge-on disks) and negligible in the opposite case, when the angle between the observer and the symmetry axes of the system is small (face-on disks). In all these spectra, the inner edge of the disk is at the ISCO radius (Novikov–Thorne model), while the outer edge is assumed at r_out = 10⁵ M. The latter is large enough that for the energy range shown in Figures 1–4 is equivalent to an infinite value. The effect of r_out is indeed to introduce a different slope of the spectrum at low energies, see, e.g., Figure 4 in Bambi & Barausse (2011a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Properties of the Backgrounds Shown in Figures 1–7: Spin Parameter a (First Column), Deformation Parameter ₃ (Second Column), Radius of the Event Horizon r_H (Third Column), Radius of the ISCO r_ISCO (Fourth Column), and Stability of the ISCO (MRS/MVS = Marginally Radially/Vertically Stable; Fifth Column)

a/M		3		rH/M		rISCO/M		ISCO
0.0		0		2		6		MRS
0.5		0		1.8660		4.2330		MRS
0.9		0		1.4359		2.3209		MRS
0.999		0		1.0447		1.1818		MRS
0.8		2		1.6 (θ = 0, π)		1.9048		MVS
0.9		1		1.4359 (θ = 0, π)		1.7499		MVS
0.9		−1		1.6176 (θ = π/2)		3.2304		MRS
1.1		−0.5		1.2801 (θ = π/2)		2.1139		MRS
1.2		−1		1.3614 (θ = π/2)		2.5036		MRS

Notes. When ₃ ≠ 0, the radius of the event horizon in general depends on the angle θ and the topology of the event horizon may be non-trivial (see Bambi & Modesto 2011). For the cases (a/M, ₃) = (0.8, 2) and (0.9, 1), there is no event horizon on the equatorial plane θ = π/2.

Download table as:  ASCIITypeset image

The Kerr nature of astrophysical BH candidates can be tested by studying the thermal spectrum of the accretion disk without the assumption of the Kerr background. The basic idea is to consider a more general spacetime, in which the compact object is characterized by a mass M, a spin parameter a, and a deformation parameter which measures possible deviation from the Kerr geometry. When the deformation parameter is set to zero, we have to recover exactly the Kerr solution. We can then compare the theoretical predictions with observational data: if the latter demand a zero deformation parameter, the Kerr BH hypothesis is verified; otherwise, we may conclude that the object is not a Kerr BH.

Let us note that this approach can be used just to check if the geometry of the spacetime is described by the Kerr solution, but it cannot really investigate the actual nature of the compact object or the exact deviations from the Kerr background. The analysis of the disk's spectrum can only probe the geometry of the spacetime till the ISCO radius; we have no information about the geometry at smaller radii and about the surface/horizon of the BH candidate. Current (and near-future) data are not so good to map the spacetime and a single deformation parameter is used to figure out if the gravitational force is stronger or weaker than the one around a Kerr BH with the same mass and spin. Actually, the typical situation is even worse and we can infer only one parameter: if we assume the Kerr background, we find the spin a, if we have also a deformation parameter, we constrain some combination of a and of the deformation parameter. However, this degeneracy can be solved with an additional measurement (Bambi 2012b, 2012c, 2012d).

When the computation of the properties of the thermal emission depends on the geometry of the background, the code calls the function metric. It is thus easy to change the spacetime. For the time being, the code uses the Johannsen–Psaltis (JP) metric, which was explicitly proposed in Johannsen & Psaltis (2011) to test the geometry around BH candidates. In Boyer–Lindquist coordinates, the metric is given by the line element (Johannsen & Psaltis 2011)

$$\begin{eqnarray} ds^2 &=& {-} \left(1 - \frac{2 M r}{\Sigma }\right) (1 + h) \, dt^2 - \frac{4 a M r \sin ^2\theta }{\Sigma } (1 + h) \, dt \, d\phi \nonumber \\ && + \frac{\Sigma (1 + h)}{\Delta + a^2 h \sin ^2\theta } \, dr^2 + \Sigma \, d\theta ^2 \nonumber\\ &&+ \left[\sin ^2\theta \left(r^2 + a^2 + \frac{2 a^2 M r \sin ^2\theta }{\Sigma } \right) \right. \nonumber\\ && +\left. \frac{a^2 (\Sigma + 2 M r) \sin ^4\theta }{\Sigma } h \right] d\phi ^2, \end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray} \Sigma &=& r^2 + a^2 \cos ^2\theta, \nonumber \\ \Delta &=& r^2 - 2 M r + a^2, \nonumber \\ h &=& \sum _{k = 0}^{\infty } \left(\epsilon _{2k} + \frac{M r}{\Sigma } \epsilon _{2k+1} \right) \left(\frac{M^2}{\Sigma }\right)^k. \end{eqnarray} \tag{ 19 }$$

This metric has an infinite number of deformation parameters _i and the Kerr solution is recovered when all the deformation parameters vanish. However, in order to reproduce the correct Newtonian limit, we have to impose ₀ = ₁ = 0, while ₂ is strongly constrained by solar system experiments. As done in other papers, it is usually sufficient to restrict the attention to the deformation parameter ₃ and set to zero all the others.

The properties of JP BHs with M, a, and ₃ have been widely discussed in the recent literature (Bambi 2011d, 2012a, 2012b, 2012c, 2012d; Bambi & Modesto 2011; Bambi et al. 2012; Johannsen & Psaltis 2012; Chen & Jing 2012; Liu et al. 2012; Krawczynski 2012). One can note that these objects may have properties fundamentally different from the ones expected for Kerr BHs. For instance, when the gravitational force on the equatorial plane turns out to be stronger/weaker than the Kerr case (this corresponds to ₃ < 0/₃ > 0), the maximum value of a/M can presumably be larger/smaller than 1, and increases/decreases as ₃ decreases/increases (Bambi 2011d). Another important feature to bear in mind is that the ISCO may not be determined by the orbital stability along the radial direction, as in Kerr, but it may be marginally vertically stable. The ISCO radius should thus be computed as described in Bambi & Barausse (2011a, 2011b). The event horizon of JP BHs has been discussed in Bambi & Modesto (2011) and the radius r_H can be found from³

$$\begin{eqnarray} &&\Delta + a^2 h \sin ^2 \theta = 0. \end{eqnarray} \tag{ 20 }$$

For non-vanishing a and h, the radius r_H is a function of the angle θ. It is also possible to show that the topology of the event horizon of these BHs may be non-trivial and that BHs with a topologically non-trivial event horizon may be created by astrophysical processes like the accretion from disk (Bambi & Modesto 2011). The fundamental properties of the JP spacetimes shown in Figures 3–7 are reported in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Examples of the observed spectrum of thin accretion disks in the JP background are shown in Figure 3, for a/M = 0.8 and ₃ = 2 (left panel) and for a/M = 1.2 and ₃ = −1 (right panel). Figure 3 shows the effect of light bending in the JP background. As in the Kerr case, the effect is more and more important as the disk's inclination angle increases and the observer approaches the equatorial plane of the system. Figure 4 shows instead the effect on the observed spectrum of the properties of the emission of the accretion disk. It compares the cases of isotropic (ϒ = 1) and limb-darkened emission (ϒ = 1/2  +  (3/4)cos ξ, ξ being the angle between the wavevector of the photon emitted by the disk and the normal to the disk surface). The left panel of Figure 4 can be compared with Figure 9 of Li et al. (2005). As noticed in Li et al. (2005), for edge-on disks the effect of limb darkening cannot be absorbed by a redefinition of $\dot{M}$ and f_col; the spectrum has indeed a more pronounced hump before the exponential cutoff. The case of non-Kerr background does not introduce any new qualitatively different feature with respect to the Kerr case already discussed in the literature.

In the case of a geometrically thick and optically thin accretion disk around a BH candidate, it should be possible to observe the so-called BH shadow; that is, a dark area over a brighter background. As the exact shape of the shadow is determined exclusively by the geometry of the spacetime, its detection can be used to investigate the nature of the BH candidate (Bambi & Freese 2009; Bambi & Yoshida 2010; Johannsen & Psaltis 2010; Bambi et al. 2012). In the case of geometrically thin and optically thick accretion disks, the image turns out to be quite different. For Schwarzschild and Kerr BHs, it has already been discussed in the literature (Luminet 1979; Fukue & Yokoyama 1988; Fukue 2003; Takahashi 2004). The generalization to non-Kerr spacetimes is straightforward (see also Krawczynski 2012). The temperature and the photon flux as detected by a distant observer are given by

$$\begin{eqnarray} T_{\rm obs} &=& g T_{\rm col}, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \mathcal {F}_{\rm obs} &=& g^4 \mathcal {F}, \end{eqnarray} \tag{ 22 }$$

where g is the redshift factor given by Equation (5). Equation (21) follows from the definition of g and from Equation (3), while Equation (22) is a consequence of Liouville's theorem I_e(ν_e)/ν³_e = I_obs(ν_obs)/ν³. Figure 5 shows the direct image of T_obs and $\mathcal {F}_{\rm obs}$ for a Schwarzcshild BH (top panels) and for a Kerr BH with a/M = 0.9 (bottom panels), in the case of an observer with inclination angle of 80°. The asymmetry with respect to the X axes is due to the effect of light bending, while that with respect to the Y axes comes from the Doppler boosting. Figure 6 shows the same images in the case when the accretion disk is in the JP spacetime with spin a/M = 0.9 and deformation parameter ₃ = −1 (top panels) and 1 (bottom panels). As we can see, the shape of the central dark area and the intensity map of the image depend clearly on the geometry of the spacetime around the compact object. Such images are definitively impossible to observe with current X-ray facilities, but they may be seen with future X-ray interferometry techniques.

The observation of the light curve when the accretion disk around a BH candidate is eclipsed by the stellar companion can also provide information about the geometry of the spacetime around the compact object (Fukue & Yokoyama 1988). While current observational facilities would already have the correct time resolution to observe this feature, so far we know only one binary system with a BH candidate showing an eclipse in the X-ray spectrum. This source is M33 X-7, which is unfortunately quite far and dim, and present detectors cannot measure its light curve with an accuracy to observe any relativistic effect (Pietsch et al. 2004).

Figure 7 shows the light curves of the ingress (left panel) and egress (right panel) during an eclipse for the four cases shown in Figures 5 and 6. The inclination angle of the observer is still i = 80°. In these simulations, I assume that the angular momentum of the companion star is parallel to that of the accretion disk, as one should probably expect from considerations on the evolution of the system. That means that the companion star occults first the approaching (blueshifted) part of the disk and then the receding (redshifted) one. The companion star is simply modeled as an obstacle with vertical edge and its atmosphere is completely neglected, so the disk's luminosity at the time t is

$$\begin{eqnarray} &&L(t) = \int h_i (X - v_{\rm c} t) I_{\rm obs}(\nu) \frac{dX dY}{D^2}, \end{eqnarray} \tag{ 23 }$$

where i = ingress or egress, v_c is the velocity of the stellar companion, and

$$\begin{eqnarray} h_{\rm ingress} (X - v_{\rm c} t) = \left\lbrace \begin{array}{@{} ll}0 & {\rm if }\; X - v_{\rm c} t \le 0, \\ \ms 1 & {\rm if }\; X - v_{\rm c} t > 0. \end{array} \right. \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} h_{\rm egress} (X - v_{\rm c} t) = \left\lbrace \begin{array}{@{} ll}1 & {\rm if }\; X - v_{\rm c} t < 0, \\ \ms 0 & {\rm if }\; X - v_{\rm c} t \ge 0. \end{array} \right. \end{eqnarray} \tag{ 25 }$$

As the shape of the light curve depends on the background metric, the possible detection of this feature can be used to constrain possible deviations from the Kerr nature of a BH candidate. However, in order to extract information on the spacetime geometry from real data, a reliable atmospheric model would be necessary (Takahashi & Watarai 2007).

If we could get accurate images of the disk, like the ones shown in Figures 5 and 6, we could determine at the same time a and ₃. However, that is surely out of reach in the near future. In the case of the light curve during an X-ray eclipse, basically we measure the slope of the curve and the asymmetry between ingress and egress. Assuming a Kerr background, we could immediately estimate a (assuming M is known independently). If we are going to test the Kerr nature of the BH candidate, we find a degeneracy between a and ₃. However, the light curve is essentially sensitive only to the effect of Doppler boosting; that is, the asymmetry of the apparent image of the disk with respect to the X = 0 axes. It is not very sensitive to the effect of light bending (responsible for the asymmetry with respect to the Y = 0 axes) or to the gravitational redshift (which produces corrections symmetric with respect to the X = 0 axes). The combination of the fit of the thermal spectrum and of the light curve during an eclipse for the same BH candidate should thus break the degeneracy between a and ₃ and allows for the identification of a limited allowed region in the spin parameter–deformation parameter plane.