1 Introduction

Astrophysical black hole candidates are dark and compact objects that can be naturally interpreted as black holes and they may be something else only in the presence of new physics. Stellar-mass black hole candidates are compact objects in X-ray binaries with a mass exceeding the maximum mass for a neutron star [1]. Supermassive black hole candidates are the huge compact bodies at the center of every normal galaxy and they turn out to be too massive, compact, and too old to be a cluster of neutron stars [2]. The non-detection of thermal radiation from the surface of these objects is also consistent with the idea that they do not have a surface but an event horizon [3, 4].

According to general relativity, the spacetime metric around black hole candidates should be well described by the Kerr solution. Initial deviations from the Kerr metric are quickly radiated away through the emission of gravitational waves [5]. The equilibrium electric charge is completely negligible for macroscopic objects [6]. The accretion disk is usually many orders of magnitude smaller than the central black hole candidate and it cannot appreciably change the geometry of the spacetime [7].

The Kerr black hole hypothesis entirely relies on the validity of general relativity and there is no clear observational confirmation that the spacetime geometry around black hole candidates is described by the Kerr solution. Moreover, general relativity has been tested only for weak gravitational fields and it is not guaranteed that its predictions still hold in the strong gravity regime.

In the past few years, there have been a significant work to study how present and future observational facilities could test black hole candidates; see e.g. Refs. [8, 9]. The most common approach is to use a method similar to the parametrized post-Newtonian (PPN) formalism [10], in which one wants to test the Schwarzschild solution in the weak field limit. In the case of black hole candidates, one employs a metric that is parametrized by a number of deformation parameters capable of describing possible deviations from the Kerr background. The deformation parameters are free quantities to be determined by observations, and a posteriori one can check whether astronomical data require vanishing deformation parameters, as it is required by the Kerr solution.

In the literature there is already a number of parametrizations suitable to test black hole candidates [11–16]. Each proposal has its own advantages and disadvantages. However, in the analysis of real data it is necessary to calculate a large number of spectra for different values of the model parameters in order to find the best fit and measure the model parameters. Typical calculations require the determination of the point of the emission on the accretion disk and of its redshift factor [17–23]. The Kerr metric has some nice properties, and eventually these calculations can be done in a reasonable time with current computational facilities. The non-Kerr metrics used to test black hole candidates do not have such nice properties and this becomes an issue when we want to measure the deformation parameters with real data.

In the present paper, we discuss a parametrization suitable for numerical calculations involving the electromagnetic spectrum of thin disks. This is our main motivation, and in particular we have in mind the continuum-fitting and the iron line methods [17–24]. We are thus interested in a metric with properties similar to the Kerr solution, and we do not look for a very general black hole metric. Our metric has the Carter constant and therefore the equations of motion are separable and of first order. More importantly, the motion along the \(\theta \)-direction is like in Kerr metric and it can be reduced to an elliptic integral, while the motion along the r-direction can be reduced to a hyper-elliptic integral. In general, this is probably the best we can have from the point of view of the accuracy and the speed of the calculations. We also note that our metric has a quite rich phenomenology. For instance, it can describe black holes with a very high Novikov–Thorne radiative efficiency, which is not the case for most (if not all) parametrizations already proposed in the literature. Such a property is quite useful when we have fast-rotating objects, which are the best candidates to test the Kerr paradigm. Another feature, which is usually absent in the other parametrization, is that our black holes can be very small. It is also worth noting that our metric has no curvature singularities outside of the event horizon.

2 Metric

As in the other parametrizations discussed in the literature [11–16], even our choice is necessarily ad hoc and it can only be motivated by our requirements, which are determined by the specific use we have in mind. The metric must clearly includes the Kerr solution as a special case. We want that there is the Carter constant, so that it is not necessary to solve the geodesic equations but the equations of motion are separable and of first order. To do this, we write the Kerr metric in Boyer–Lindquist coordinates and we promote the constant M to some functions \(m_i (r)\) which depend on the radial coordinate only. The line element reads

$$\begin{aligned} \text {d}s^2= & {} - \left( 1 - \frac{2 m_1 r}{\Sigma }\right) \text {d}t^2 - \frac{4 a m_1 r \sin ^2\theta }{\Sigma } \text {d}t\text {d}\phi \nonumber \\&+ \frac{\Sigma }{\Delta _2} \text {d}r^2 + \Sigma \text {d}\theta ^2 \nonumber \\&+ \left( r^2 + a^2 + \frac{2 a^2 m_1 r \sin ^2\theta }{\Sigma }\right) \sin ^2\theta \text {d}\phi ^2 , \end{aligned}$$
(1)

where \(\Sigma = r^2 + a^2 \cos ^2\theta \), \(\Delta _2 = r^2 - 2 m_2 r + a^2\), \(m_1 = m_1(r)\), and \(m_2 = m_2(r)\). In \(g_{tt}\), \(g_{t\phi }\), and \(g_{\phi \phi }\), M has been replaced by the same function \(m_1\), because otherwise we lose the separability of the equations of motion (see the next section). The mass M in \(g_{rr}\) has been replaced by the function \(m_2\), which (in general) may be different by \(m_1\) without affecting our requirement.

In the metric in Eq. (1), the radius of the event horizon, \(R_\mathrm{H}\), is given by the largest root in \(g^{rr} = 0\), namely \(\Delta _2 = 0\), where only \(m_2\) is involved, not \(m_1\). On the other hand, the Killing horizon is given by the largest root in

$$\begin{aligned} g_{tt} g_{\phi \phi } - g_{t\phi }^2 = 0 . \end{aligned}$$
(2)

If \(m_1 = m_2\), an extension of the rigidity theorem holds and event horizon and Killing horizon coincide. In the general case, with \(m_1 \ne m_2\), this may not be true.

We note that the metric in Eq. (1) reduces to the Kerr–Newman solution when

$$\begin{aligned} m_1 = m_2 = M - \frac{Q^2}{2 r} , \end{aligned}$$
(3)

where Q is the electric charge of the black hole. Our metric can also describe a large class of quantum gravity inspired black hole solutions [25–29] and some regular black holes [30, 31]. For instance, in the non-commutative inspired black holes of Refs. [25–28], one has

$$\begin{aligned} m_1(r) = m_2(r) = \frac{\gamma (3/2;r^2/4l_0^2)}{\Gamma (3/2)} \, M, \end{aligned}$$
(4)

where \(\gamma (3/2;r^2/4l_0^2)\) is the lower incomplete Gamma function, \(\Gamma (3/2) = \sqrt{\pi }/2\) is the Gamma function at 3 / 2, and \(l_0\) is the non-commutativity length scale of the theory. In the weakly non-local theories of gravity of Ref. [29], the black hole solutions have

$$\begin{aligned}&m_1(r) = m_2(r) \nonumber \\&\quad = \frac{2M}{\pi } \int _0^r \text {d}x \, x^2 \int ^\infty _0 \text {d}k \, k^2 \frac{\sin kr}{kr} V\left( -\frac{k^2}{\Lambda ^2}\right) , \end{aligned}$$
(5)

where V is the model-dependent form factor and \(\Lambda \) is the scale of the theory. The Bardeen metric [30, 31] has

$$\begin{aligned} m_1(r) = m_2(r) = \frac{r^3}{\left( r^2 + g^2\right) ^{3/2}} \, M . \end{aligned}$$
(6)

If the Bardeen solution is derived from Einstein gravity coupled to a non-linear electrodynamics field [32], g is the magnetic charge of the black hole. In all these examples, we always have \(m_1 = m_2\), but in the following we will consider the general case without this condition when it is not required the exact expression of \(m_1\) and \(m_2\).

Let us write \(m_1\) and \(m_2\) in the following form:

$$\begin{aligned} m_i = M \sum _{k=0}^{\infty } a_{ik} \left( \frac{M}{r}\right) ^k . \end{aligned}$$
(7)

In the weak field regime, \(M/r \ll 1\) and can be used as an expansion parameter. The metric coefficient \(g_{tt}\) and \(g_{rr}\) become

$$\begin{aligned} g_{tt}= & {} - \left[ 1 - a_{10} \frac{2 M}{r} - a_{11} \frac{2 M^2}{r^2} + \cdots \right] , \end{aligned}$$
(8)
$$\begin{aligned} g_{rr}= & {} 1 + a_{20} \frac{2 M}{r} + \cdots . \end{aligned}$$
(9)

When cast in Schwarzschild coordinates, the PPN metric reduces to

$$\begin{aligned} g_{tt}= & {} - \left[ 1 - \frac{2 M}{r} + \left( \beta _\mathrm{PPN} - \gamma _\mathrm{PPN} \right) \frac{2 M^2}{r^2} + \cdots \right] , \end{aligned}$$
(10)
$$\begin{aligned} g_{rr}= & {} 1 + \gamma \frac{2 M}{r} + \cdots , \end{aligned}$$
(11)

and Solar System experiments require that \(\beta _\mathrm{PPN}\) and \(\gamma _\mathrm{PPN}\) are 1 with an accuracy at the level of \(10^{-5} - 10^{-4}\) [33]. Within our parametrization and \(m_1\) and \(m_2\) given by Eq. (7), we can always choose \(a_{10} = 1\) (if \(a_{10} \ne 1\), we just redefine M). \(a_{11} \approx 0\) and \(a_{20} \approx 1\) are constrained by Solar System experiments. The first unconstrained coefficients are thus \(a_{12}\) and \(a_{21}\).

In the next sections, we will focus the attention on the following choice of \(m_1\) and \(m_2\):

$$\begin{aligned} m_1 = m_2 = M \left( 1 + \alpha \frac{M^2}{r^2} + \beta \frac{M^3}{r^3} \right) , \end{aligned}$$
(12)

where \(\alpha \) and \(\beta \) are the two deformation parameters of our metric. This choice is obtained from the truncation of the general expression in Eq. (7); that is, we consider the two leading order terms without Solar System constraints and we neglect higher order corrections. The Kerr solution is recovered when \(\alpha = \beta = 0\). With the choice in Eq. (12), there are no naked singularities in the region outside of the black hole. In the appendix, we report the expressions of the invariants R, \(R^{\mu \nu } R_{\mu \nu }\), and \(R^{\mu \nu \rho \sigma } R_{\mu \nu \rho \sigma }\). These invariants are everywhere regular except at \(r=0\).

Our final goal is to have a parametrization suitable for the numerical calculations necessary in the study of the electromagnetic spectrum of thin disks. The background metric enters the following calculations:

  1. 1.

    The motion of the particle in the disk. In the Novikov–Thorne model [34], the disk is on the equatorial plane orthogonal to the black hole spin, and the particles of the accretion disk follow nearly geodesic circular orbits on the equatorial plane.

  2. 2.

    The photon trajectories from the emission point in the disk to the detection point at infinity. Actually, it is not strictly necessary to compute the exact photon trajectories, but it is necessary to connect the emission point on the disk to that on the image plane of the distant observer [35, 36].

3 Motion of massive particles in the equatorial plane

In the calculations of the motion of the particles in the disk, only the metric coefficient \(g_{tt}\), \(g_{t\phi }\), and \(g_{\phi \phi }\) are involved. So we need to specify \(m_1\), which will be assumed to have the form in Eq. (12), while \(m_2\) is actually irrelevant in this part.

The angular velocity of equatorial circular orbits is [17, 37]

$$\begin{aligned} \Omega _\pm= & {} \frac{- \left( \partial _r g_{t\phi }\right) \pm \sqrt{\left( \partial _r g_{t\phi }\right) ^2 - \left( \partial _r g_{tt}\right) \left( \partial _r g_{\phi \phi }\right) }}{\partial _r g_{\phi \phi }} \nonumber \\= & {} \frac{\tilde{m}_1^{1/2}}{r^{3/2} \pm a \tilde{m}_1^{1/2}} , \end{aligned}$$
(13)

where here and in the next formulas the upper sign refer to corotating orbits, while the lower sign to counterrotating orbits. Moreover, we have introduced \(\tilde{m}_1\), which is defined by

$$\begin{aligned} \tilde{m}_1 = M \left( 1 + 3 \alpha \frac{M^2}{r^2} + 4 \beta \frac{M^3}{r^3} \right) \end{aligned}$$
(14)

When \(\alpha = \beta = 0\), we clearly recover the well-known Kerr result [37]

$$\begin{aligned} \Omega _\pm ^\mathrm{Kerr} = \frac{M^{1/2}}{r^{3/2} \pm a M^{1/2}} . \end{aligned}$$
(15)

Following the standard prescription to find the equatorial circular orbits [17, 37], the specific energy and specific angular momentum of the particles of the gas are

$$\begin{aligned} E= & {} \frac{r^{3/2} - 2 m_1 r^{1/2} \pm a \tilde{m}_1^{1/2}}{r^{3/4} \sqrt{r^{3/2} - 2 m_1 r^{1/2} - \tilde{m}_1 r^{1/2} \pm 2 a \tilde{m}_1^{1/2}}} , \end{aligned}$$
(16)
$$\begin{aligned} L_z= & {} \pm \frac{\tilde{m}_1^{1/2} \left( r^2 \mp 2 a m_1^{1/2} \tilde{m}_1^{-1/2} r^{1/2} + a^2\right) }{r^{3/4} \sqrt{r^{3/2} - 2 m_1 r^{1/2} - \tilde{m}_1 r^{1/2} \pm 2 a \tilde{m}_1^{1/2}}} . \end{aligned}$$
(17)

We recover the correct Kerr limit for \(\alpha = \beta = 0\)

$$\begin{aligned} E= & {} \frac{r^{3/2} - 2 M r^{1/2} \pm a M^{1/2}}{r^{3/4} \sqrt{r^{3/2} - 3 M r^{1/2} \pm 2 a M^{1/2}}} , \end{aligned}$$
(18)
$$\begin{aligned} L_z= & {} \pm \frac{M^{1/2} \left( r^2 \mp 2 a M^{1/2} r^{1/2} + a^2\right) }{r^{3/4} \sqrt{r^{3/2} - 3 M r^{1/2} \pm 2 a M^{1/2}}}. \end{aligned}$$
(19)

When

$$\begin{aligned} r^{3/2} - 2 m_1 r^{1/2} - \tilde{m}_1 r^{1/2} \pm 2 a \tilde{m}_1^{1/2} = 0 , \end{aligned}$$
(20)

the denominator in Eqs. (17) and (16) vanishes and the particle has infinite energy. Equation (20) defines the radius of the photon orbit, \(R_\mathrm{photon}\), which is the minimum radius for circular orbits (at smaller radii, there are no circular orbits).

In the spectrum of thin disks, the inner edge of the disk plays a crucial role and it is eventually the true key-ingredient in both the continuum-fitting and the iron line methods. In the Novikov–Thorne model, the inner edge of the disk is at the innermost stable circular orbit (ISCO). In general, one has to check the orbital stability along both the radial and the vertical directions [17, 38]. In our spacetime we have checked that the ISCO radius is only determined by the orbital stability along the radial direction (like in the Kerr metric), at least for not too large values of the deformation parameters. In this case, the ISCO radius corresponds to the minimum of the specific energy

$$\begin{aligned} \frac{\text {d}E}{\text {d}r} = 0 \Rightarrow r = R_\mathrm{ISCO}. \end{aligned}$$
(21)

Unfortunately, it seems there is not a compact analytic expression for \(R_\mathrm{ISCO}\) as in the Kerr metric.

The so-called Novikov–Thorne radiative efficiency, which is the actual quantities measured by the continuum-fitting method [20], is

$$\begin{aligned} \eta _\mathrm{NT} = 1 - E_\mathrm{ISCO} , \end{aligned}$$
(22)

where \(E_\mathrm{ISCO}\) is the specific energy of a test particle at the ISCO radius, namely Eq. (17) evaluated at \(r = R_\mathrm{ISCO}\).

The 4-velocity of a particle in the disk is given by \(u^\mu _\mathrm{e} = (u^t_\mathrm{e}, 0, 0, u^\phi _\mathrm{e})\), where \(u^\phi _\mathrm{e} = \Omega u^t_\mathrm{e}\) by definition of \(\Omega \). From the normalization condition \(g_{\mu \nu } u^\mu _\mathrm{e} u^\nu _\mathrm{e} = -1\), we get the expression for \(u^t_\mathrm{e}\)

$$\begin{aligned} u^t_\mathrm{e}= & {} \frac{1}{\sqrt{-g_{tt} - 2g_{t\phi }\Omega - g_{\phi \phi }\Omega ^2}} \nonumber \\= & {} \frac{r^{3/2}_\mathrm{e} + a \tilde{m}_1^{1/2}}{r^{1/2}_\mathrm{e} \sqrt{r^2_\mathrm{e} - 2 m_1 r_\mathrm{e} - \tilde{m}_1 r_\mathrm{e} + 2 a \tilde{m}_1^{1/2} r^{1/2}_\mathrm{e}}} , \end{aligned}$$
(23)

which reduces to the correct Kerr case for \(\alpha = \beta = 0\)

$$\begin{aligned} u^t_\mathrm{e} = \frac{r^{3/2}_\mathrm{e} + a M^{1/2}}{r^{1/2}_\mathrm{e} \sqrt{r^2_\mathrm{e} - 3 M r_\mathrm{e} + 2 a M^{1/2} r^{1/2}_\mathrm{e}}}. \end{aligned}$$
(24)

The 4-velocity of the particles in the disk is necessary to calculate another important quantity in the calculation of the spectrum of thin disks, namely the redshift factor g

$$\begin{aligned} g= \frac{\nu _\mathrm{o}}{\nu _\mathrm{e}} = \frac{u^\mu _\mathrm{o} k_\mu }{u^\nu _\mathrm{e} k_\nu } = \frac{\sqrt{-g_{tt} - 2g_{t\phi }\Omega - g_{\phi \phi }\Omega ^2}}{1 - \xi \Omega } , \end{aligned}$$
(25)

where \(\nu _\mathrm{o}\) and \(\nu _\mathrm{e}\) are, respectively, the photon frequency as measured by the distant observer and the emitter, \(u^\mu _\mathrm{o} = (1,0,0,0)\) is the 4-velocity of the observer, \(k^\mu \) is the 4-momentum of the photon, and \(\xi = - k_\phi /k_t\) is a constant of motion along the photon trajectory (as a consequence of the fact the spacetime is stationary and axisymmetric). Within our parametrization, we find (for \(a \ge 0\))

$$\begin{aligned} g = \frac{r^{1/2}_\mathrm{e} \sqrt{r^2_\mathrm{e} - 2 m_1 r_\mathrm{e} - \tilde{m}_1 r_\mathrm{e} + 2 a \tilde{m}_1^{1/2} r^{1/2}_\mathrm{e}}}{r^{3/2}_\mathrm{e} + a \tilde{m}_1^{1/2} - \tilde{m}_1^{1/2} \xi } . \end{aligned}$$
(26)

The correct Kerr limit is again recovered for \(\alpha = \beta = 0\)

$$\begin{aligned} g = {\frac{r_\mathrm{e}^{1/2} \sqrt{r_\mathrm{e}^2 - 3 M r_\mathrm{e} + 2 a M^{1/2} r^{1/2}_\mathrm{e}}}{r^{3/2}_\mathrm{e} + a M^{1/2} - M^{1/2} \xi }} . \end{aligned}$$
(27)

In the end, we have \(g = g(r_\mathrm{e},\xi )\). Since \(r_\mathrm{e} = r_\mathrm{e}(\xi ,q)\), where \(q^2={\mathcal Q}/E^2\) and \({\mathcal Q}\) is the Carter constant (see Sect. 5), eventually we have [35]

$$\begin{aligned} g = g (\xi , q) . \end{aligned}$$
(28)

Figures 1, 2, 3, 4, and 5 show the contour levels of the radius of the event horizon \(R_\mathrm{H}\) (top left panels), of the photon radius \(R_\mathrm{photon}\) (top right panels), of the ISCO radius \(R_\mathrm{ISCO}\) (bottom left panels), and of the Novikov–Thorne radiative efficiency \(\eta _\mathrm{NT}\) (bottom right panels) for \(\beta = 0\) (Fig. 1), \(-0.2\) (Fig. 2), \(-0.5\) (Fig. 3), 0.2 (Fig. 4), and 0.5 (Fig. 5).

The light-green regions are the parameter space with black holes, the white regions are those with naked singularities (no real solution of the equation \(r^2 - 2 m_2 r + a^2 = 0\)). \(\alpha , \beta > 0\) make the gravitational force at small radii stronger (they “increase” the value of the effective mass), while \(\alpha , \beta < 0\) make it weaker. The result is that for \(\beta = 0\), \(-0.2\), and \(-0.5\), there may be a naked singularity for \(\alpha < 0\) because in those cases the gravitational force is not strong enough to create an event horizon. For \(\beta = 0.5\), there are no naked singularities in the plots because at small radii the dominant term is \(\beta M^3/r^3\) and it is positive (namely gravity is strong and there is an event horizon).

It is worth noting that the Novikov–Thorne radiative efficiency \(\eta _\mathrm{NT}\) of our black holes may exceed the maximum value for a Kerr black hole \(\eta _\mathrm{NT}^\mathrm{max} \approx 0.42\). This is not the case for most parametrizations proposed in the literature. If a black hole candidate has a high radiative efficiency, within the parametrizations in the literature it is quite automatic that the deviations from Kerr must be small, or otherwise it is impossible to reproduce the spectrum. In our case, this is not true, which means that our parametrization includes deviations from the Kerr solutions that are usually not taken into account. In a similar way, our black holes can be very small (the radius of the event horizon is small) and may also appear small (the photon capture radius is small). This is not the case in the other parametrizations.

Fig. 1
figure 1

Contour levels of the radius of the event horizon \(R_\mathrm{H}\) (top left panel), of the photon radius \(R_\mathrm{photon}\) (top right panel), of the ISCO radius \(R_\mathrm{ISCO}\) (bottom left panel), and of the Novikov–Thorne radiative efficiency \(\eta _\mathrm{NT}\) (bottom right panel) for \(\beta = 0\). \(m_1\) and \(m_2\) are given by Eq. (12). The spacetimes in the white region have no black hole but a naked singularity and they have not been studied. \(a_* = a/M\) is the dimensionless spin parameter. \(R_\mathrm{H}\), \(R_\mathrm{photon}\), and \(R_\mathrm{ISCO}\) in units with \(M=1\). See the text for more details

Full size image
Fig. 2
figure 2

As in Fig. 1 for \(\beta = -0.2\)

Full size image
Fig. 3
figure 3

As in Fig. 1 for \(\beta = -0.5\)

Full size image
Fig. 4
figure 4

As in Fig. 1 for \(\beta = 0.2\)

Full size image
Fig. 5
figure 5

As in Fig. 1 for \(\beta = 0.5\)

Full size image

4 Hamilton–Jacobi equation

The second ingredient necessary in the calculation of the spectrum of a thin disk is the determination of the photon trajectories from the point of emission to the point of detection. In a general spacetime, this is done by solving the geodesic equations, which are second order partial differential equations in the coordinates of the spacetime. The Kerr metric has the Carter constant and the equations of motion are separable and of first order. More importantly, the motion along the \(\theta \) and r directions can be reduced to elliptic integrals. The result is that numerical calculations can be faster and more accurate. In any non-trivial extension of the Kerr metric, this is not possible. However, we can have a metric in which the motion along the \(\theta \) and r directions can be reduced to hyper-elliptic integrals, with similar advantages of the Kerr solution.

The starting point is the Hamilton–Jacobi equation

$$\begin{aligned} 2 \frac{\partial S}{\partial \tau } = g^{\mu \nu } \frac{\partial S}{\partial x^\mu } \frac{\partial S}{\partial x^\nu } . \end{aligned}$$
(29)

Assuming in this section the more general case with \(m_1\) and \(m_2\) not necessarily the same, \(g^{\mu \nu }\) is

$$\begin{aligned} \left( \frac{\partial }{\partial s}\right) ^2= & {} - \frac{A_1}{\Sigma \Delta _1} \left( \frac{\partial }{\partial t}\right) ^2 - \frac{4 a m_1 r}{\Sigma \Delta _1} \left( \frac{\partial }{\partial t}\right) \left( \frac{\partial }{\partial \phi }\right) \nonumber \\&+ \frac{\Delta _2}{\Sigma } \left( \frac{\partial }{\partial r}\right) ^2 + \frac{1}{\Sigma } \left( \frac{\partial }{\partial \theta }\right) ^2 \nonumber \\&+ \frac{\Delta _1 - a^2 \sin ^2\theta }{\Sigma \Delta _1 \sin ^2\theta } \left( \frac{\partial }{\partial \phi }\right) ^2 , \end{aligned}$$
(30)

where \(\Delta _1 = r^2 - 2 m_1 r + a^2\) and \(A_1 = \left( r^2 + a^2\right) ^2 - a^2 \Delta _1 \sin ^2\theta \).

We can then proceed as in the Kerr case, looking for a solution of the Hamilton–Jacobi equation of the form (see Ref. [39], Chapter 7, Section 62 for the details)

$$\begin{aligned} S= & {} - \frac{1}{2} \delta \tau - E t + L_z \phi S_r(r) + S_\theta (\theta ) . \end{aligned}$$
(31)

The solution for S is

$$\begin{aligned} S= & {} - \frac{1}{2} \delta \tau - E t + L_z \phi \nonumber \\&+ \int ^r \pm dr' \sqrt{\frac{R(r')}{\Delta _1 \Delta _2}} + \int ^\theta \pm d\theta ' \sqrt{\Theta (\theta ')} , \end{aligned}$$
(32)

where \(\delta = 1\) \((\delta = 0)\) for time-like (null) geodesics, R(r) and \(\Theta (\theta )\) are given by

$$\begin{aligned} R(r)= & {} \left[ \left( r^2 + a^2\right) E - a L_z\right] ^2 \nonumber \\&- \Delta _1 \left[ {\mathcal Q} + \left( L_z - a E\right) ^2 + \delta r^2\right] , \end{aligned}$$
(33)
$$\begin{aligned} \Theta (\theta ) = {\mathcal Q} - \left[ a^2 \left( \delta - E^2\right) + L_z^2 \csc ^2\theta \right] \cos ^2\theta , \end{aligned}$$
(34)

and the signs \(\pm \) in (32) depend on the photon direction and they change at the turning points [39]. \({\mathcal Q}\) is the Carter constant, which reduces to the Carter constant of the Kerr metric when \(\alpha =\beta =0\).

The equations of motion can be obtained by setting to zero the partial derivatives of S with respect to the four constants of motion, \(\delta \), E, \(L_z\), and \({\mathcal Q}\). From \(\partial S/\partial {\mathcal Q}=0\) we get

$$\begin{aligned} \int ^r \pm \text {d}r' \sqrt{\frac{\Delta _1}{\Delta _2 R}}= & {} \int ^\theta \pm {\frac{\text {d}\theta '}{\sqrt{\Theta }}} . \end{aligned}$$
(35)

From \(\partial S/\partial \delta =0\), \(\partial S/\partial E=0\), and \(\partial S/\partial L_z=0\), we find, respectively,

$$\begin{aligned} \tau= & {} \int ^r \text {d}r' r'^2 \sqrt{\frac{\Delta _1}{\Delta _2 R}} + a^2 \int ^\theta \text {d}\theta ' \frac{\cos ^2\theta '}{\sqrt{\Theta }} , \end{aligned}$$
(36)
$$\begin{aligned} t= & {} \tau E \nonumber \\&+ 2 \int ^r \frac{m_1(r') \text {d}r'}{\sqrt{\Delta _1 \Delta _2 R}} \left[ r'^2 E - a \left( L_z - a E\right) \right] r' , \end{aligned}$$
(37)
$$\begin{aligned} \phi= & {} a \int ^r \frac{\text {d}r'}{\sqrt{\Delta _1 \Delta _2 R}} \left[ \left( r'^2 + a^2\right) E - a L_z\right] \nonumber \\&+ \int ^\theta \frac{\text {d}\theta '}{\sqrt{\Theta }} \left( L_z \csc ^2\theta ' - a E\right) . \end{aligned}$$
(38)

As in the Kerr metric [39], it is straightforward to verify that the system of Eqs. (35)–(38) is equivalent to the set of equations

$$\begin{aligned} \Sigma ^2 \dot{r}^2= & {} \frac{\Delta _2}{\Delta _1} R , \end{aligned}$$
(39)
$$\begin{aligned} \Sigma ^2 \dot{\theta }^2= & {} \Theta , \end{aligned}$$
(40)
$$\begin{aligned} \Sigma \dot{\phi }= & {} \frac{1}{\Delta _1} \left[ 2 a m_1 r E + \left( \Sigma - 2 m_1 r\right) L_z \csc ^2\theta \right] , \end{aligned}$$
(41)
$$\begin{aligned} \Sigma \dot{t}= & {} \frac{1}{\Delta _1} \left( A_1 E - 2 a m_1 r L_z\right) . \end{aligned}$$
(42)

5 Motion of massless particles from the disk to the observer

We want now to study the motion of the photons from the point of emission in the disk to the point of detection at infinity. We assume the choice in Eq. (12) and we define \(m = m_1 = m_2\). For null geodesics, \(\delta =0\) and we use the parameters \(\xi = -k_\phi /k_t\) and \(q^2 = {\mathcal Q}/k_t^2\), where \(k_t=-E\) and \(k_\phi = L_z\). This choice is useful for \(\delta =0\) because the photon trajectories are independent of the photon energy. We have

$$\begin{aligned} \tilde{R}= & {} \frac{R}{E^2} = r^4 + \left( a^2 - \xi ^2 - q^2\right) r^2 \nonumber \\&+ 2 m \left[ q^2 + \left( \xi - a\right) ^2\right] r - a^2 q^2 , \end{aligned}$$
(43)
$$\begin{aligned} \tilde{\Theta }= & {} \frac{\Theta }{E^2} = q^2 + a^2 \cos ^2\theta - \xi ^2 \cot ^2\theta . \end{aligned}$$
(44)

The relation between the parameters \((\xi ,q)\) and the celestial coordinates (X, Y) of the image as seen by an observer at infinity is the same as in Kerr (because it requires \(r\rightarrow \infty \)) [39]

$$\begin{aligned} X= & {} \lim _{r \rightarrow \infty } \left( \frac{r p^{(\phi )}}{p^{(t)}}\right) = \xi \csc \theta _\mathrm{o} , \end{aligned}$$
(45)
$$\begin{aligned} Y= & {} \lim _{r \rightarrow \infty } \left( \frac{r p^{(\theta )}}{p^{(t)}}\right) \nonumber \\= & {} \pm \sqrt{q^2 + a^2 \cos ^2\theta _\mathrm{o} - \xi ^2 \cot ^2\theta _\mathrm{o}} , \end{aligned}$$
(46)

where \(\theta _\mathrm{o}\) is the angular coordinate of the observer at infinity.

In the calculation of the photons from the disk to the distant observer, we are only interested in the motion on the \((r,\theta )\)-plane. The master equation is

$$\begin{aligned} \int ^{r_\mathrm{o}}_{r_\mathrm{e}} \frac{\text {d}r'}{\sqrt{\tilde{R}}} = \int ^{\theta _\mathrm{o}}_{\pi /2} \frac{\text {d}\theta '}{\sqrt{\tilde{\Theta }}} , \end{aligned}$$
(47)

where \(r_\mathrm{o}\) is the radial coordinate of the observer at infinity, \(r_\mathrm{e}\) is the radius of the emission point on the disk, and \(\theta _\mathrm{e} = \pi /2\) because the disk is in the equatorial plane.

Since \(\tilde{\Theta }\) is independent of \(\alpha \) and \(\beta \), one can use the same calculation technique as in Kerr [35, 39–41]. The integral can be transformed to

$$\begin{aligned} \int ^{\theta _\mathrm{o}}_{\pi /2} \frac{\text {d}\theta '}{\sqrt{\tilde{\Theta }}} = C_\theta F[\psi _\theta (\pi /2), \kappa _\theta ] , \end{aligned}$$
(48)

where F is an elliptic integral of the first kind with argument \(\psi _\theta \) and modulus \(\kappa _\theta \), while \(C_\theta \), \(\psi _\theta \), and \(\kappa _\theta \) are functions of \(\xi \) and q.

In the Kerr case, even the integral in r can be reduced to an elliptic integral of the first kind. This is not possible in a non-trivial generalization of the Kerr metric, because in the Kerr metric the function \(\tilde{R}\) in Eq. (47) is already a polynomial of fourth order. However, the integral can be transformed to a hyper-elliptic integral; see Ref. [42]. The calculations are somewhat more difficult, but the procedure is well known. In the case of our default choice (12), we have

$$\begin{aligned} \int ^{r_\mathrm{o}}_{r_\mathrm{e}} \frac{\text {d}r'}{\sqrt{\tilde{R}}}= & {} \int ^{r_\mathrm{o}}_{r_\mathrm{e}} \frac{r' \, \text {d}r'}{\sqrt{P_6(r')}} \end{aligned}$$
(49)

where \(P_6(r)\) is a polynomial of order 6,

$$\begin{aligned} P_6(r)= & {} r^6 + \left( a^2 - \xi ^2 - \eta \right) r^4 + 2 M \left[ q^2 + \left( \xi - a\right) ^2\right] r^3 \nonumber \\&- a^2 q^2 r^2 + 2 \alpha M^3 \left[ q^2 + \left( \xi - a\right) ^2\right] r \nonumber \\&+ 2 \beta M^4 \left[ q^2 + \left( \xi - a\right) ^2\right] . \end{aligned}$$
(50)

In the end, it is possible to write Eq. (49) as a function of \(\xi \), q, and \(r_\mathrm{e}\) and solve Eq. (47) in terms of \(r_\mathrm{e}\):

$$\begin{aligned} r_\mathrm{e} = r_\mathrm{e}(\xi ,q) \end{aligned}$$
(51)
Fig. 6
figure 6

Impact of the deformation parameter \(\alpha \) on the thermal spectrum of thin disks. Top left panel spectra in Kerr spacetimes with different values of the spin parameter. Top right panel spectra in spacetimes with spin parameter \(a_* = 0.7\) and different values of the deformation parameter \(\alpha \) (\(\beta = 0\)). Bottom panels spectra in non-Kerr spacetimes with \(\alpha = 0.2\) and different values of the spin parameter \(a_*\). The values of the other model parameters are: mass \(M = 10\) \(M_\odot \), mass accretion rate \(\dot{M} = 2 \times 10^{18}\) g s\(^{-1}\), distance \(D = 10\) kpc, viewing angle \(i = 45^\circ \), color factor \(f_\mathrm{col}=1.6\), and \(\Upsilon = 1\). Flux density \(N_{E_\mathrm{obs}}\) in photons keV\(^{-1}\) cm\(^{-2}\) s\(^{-1}\), and photon energy \(E_\mathrm{obs}\) in keV. See the text for more details

Full size image

6 Transfer function

The calculation of the spectrum of thin disks can be conveniently split into two parts: one concerning the background metric, and another part related to the astrophysical model. In this way, we can focus our attention on the relativistic effects determined by the background metric, which can be later combined with the astrophysical models already discussed in the literature. This can be done by introducing the transfer function f, which takes into account all the relativistic effects (gravitational redshift, Doppler boosting, light bending) [35, 36].

The observed flux is

$$\begin{aligned} F_\mathrm{o} (\nu _\mathrm{o}) = \int I_\mathrm{o}(\nu _\mathrm{o}) \text {d}\tilde{\Omega } = \int g^3 I_\mathrm{e}(\nu _\mathrm{e}) \text {d}\tilde{\Omega } , \end{aligned}$$
(52)

where \(I_\mathrm{o}\) and \(I_\mathrm{e}\) are, respectively, the specific intensities of the radiation detected by the distant observer and the specific intensities of the radiation as measured by the emitter, \(\text {d}\tilde{\Omega } = \text {d}X \text {d}Y/r^2_\mathrm{o}\) is the element of the solid angle subtended by the image of the disk on the observer’s sky, \(r_\mathrm{o}\) is the distance of the observer from the source, and \(I_\mathrm{o} = g^3 I_\mathrm{e}\) follows from Liouville’s theorem. \(I_\mathrm{e}\) is a function of the emitted frequency \(\nu _\mathrm{e}\), but also of the emission radius \(r_\mathrm{e}\) and of the direction of the photon emission, which can be described by the polar angle \(n_\mathrm{e}\) of the emitted photon with respect to the normal of the disk in the rest frame of the gas.

It is convenient to introduce the relative redshift \(g^* = g^* (r_\mathrm{e},\theta _\mathrm{o})\), defined by

$$\begin{aligned} g^* = \frac{g - g_\mathrm{min}}{g_\mathrm{max} - g_\mathrm{min}} , \end{aligned}$$
(53)

which ranges from 0 to 1. Here \(g_\mathrm{max}=g_\mathrm{max}(r_\mathrm{e},\theta _\mathrm{o})\) and \(g_\mathrm{min}=g_\mathrm{min}(r_\mathrm{e},\theta _\mathrm{o})\) are, respectively, the maximum and the minimum values of g for the photons emitted from the radial coordinate \(r_\mathrm{e}\) and detected by a distant observer with polar coordinate \(\theta _\mathrm{o}\). The observed flux can now be rewritten as

$$\begin{aligned} F_\mathrm{o} (\nu _\mathrm{o})= & {} \frac{1}{r^2_\mathrm{o}} \int _{R_\mathrm{ISCO}}^{\infty } \int _0^1 \pi r_\mathrm{e} \frac{ g^2}{\sqrt{g^* (1 - g^*)}} \nonumber \\&\times f(g^*,r_\mathrm{e},\theta _\mathrm{o}) I_\mathrm{e}(\nu _\mathrm{e},r_\mathrm{e},n_\mathrm{e}) \, \text {d}g^* \, \text {d}r_\mathrm{e} \end{aligned}$$
(54)

where f is the transfer function and takes into account all the relativistic effects determined by the background metric

$$\begin{aligned} f(g^*,r_\mathrm{e},\theta _\mathrm{o}) = \frac{1}{\pi r_\mathrm{e}} g \sqrt{g^* (1 - g^*)} \left| \frac{\partial \left( X,Y\right) }{\partial \left( g^*,r_\mathrm{e}\right) } \right| . \end{aligned}$$
(55)

Since

$$\begin{aligned} \left| \frac{\partial \left( X,Y\right) }{\partial \left( g^*,r_\mathrm{e}\right) } \right| = \frac{q \left( g_\mathrm{max} - g_\mathrm{min}\right) }{Y \sin \theta _\mathrm{o}} \left| \frac{\partial \left( \xi ,q\right) }{\partial \left( g,r_\mathrm{e}\right) } \right| , \end{aligned}$$
(56)

the calculation of the transfer function f requires the evaluation of the Jacobian

$$\begin{aligned} \left| \frac{\partial \left( \xi ,q\right) }{\partial \left( g,r_\mathrm{e}\right) } \right| = \left| \frac{\partial \xi }{\partial g} \frac{\partial q}{\partial r_\mathrm{e}} - \frac{\partial q}{\partial g} \frac{\partial \xi }{\partial r_\mathrm{e}}\right| , \end{aligned}$$
(57)

which can be done numerically from Eqs. (28) and (51) in the same way as in the Kerr case [36].

Lastly, if the transfer function depends on the emission angle \(n_\mathrm{e}\), this is given by [36]

$$\begin{aligned} \cos (n_\mathrm{e}) = - \frac{n^\mu k_\mu }{u^\nu _\mathrm{e} k_\nu } = \frac{q g}{r_\mathrm{e}} , \end{aligned}$$
(58)

because \(n^\mu = \left( 0, 0, 1/r_\mathrm{e}, 0\right) \), \(u^\nu _\mathrm{e} k_\nu = k_t/g\), and \(q = - k_\theta /k_t\).

Fig. 7
figure 7

Thermal spectrum of thin disks in non-Kerr spacetimes with \(\alpha = -0.3\) and \(\beta = 0.2\). The other model parameters are the same as in Fig. 6. These spectra can be understood in terms of the Novikov–Thorne radiative efficiency. See the text for more details

Full size image

7 Thermal spectrum of thin disks

While it is not the purpose of this paper to study the observational implications of our metric and to constrain \(\alpha \) and \(\beta \), it is useful to understand the impact of these deformation parameters on the spectrum of a black hole. To do this, we consider the thermal spectrum of a thin disk. In this case, the specific intensities of the radiation in the rest frame of the gas is (for more details, see e.g. [18] and the references therein)

$$\begin{aligned} I_\mathrm{e}(\nu _\mathrm{e},r_\mathrm{e},n_\mathrm{e}) = \frac{2 h \nu ^3_\mathrm{e}}{c^2} \frac{1}{f_\mathrm{col}^4} \frac{\Upsilon (n_\mathrm{e})}{\exp \left( \frac{h \nu _\mathrm{e}}{k_\mathrm{B} T_\mathrm{col}(r)}\right) - 1} , \end{aligned}$$
(59)

where h is the Planck constant, c is the speed of light, \(k_\mathrm{B}\) is the Boltzmann constant, \(\Upsilon (n_\mathrm{e})\) is a function that depends on the emission model (for example, \(\Upsilon = 1\) for isotropic emission and \(\Upsilon = \frac{1}{2} + \frac{3}{4} n_\mathrm{e}\) for limb-darkened emission), and \(f_\mathrm{col} \approx 1.6\) is the color (or hardening) factor. The color temperature is \(T_\mathrm{col}(r) = f_\mathrm{col} T_\mathrm{eff} (r)\), where \(T_\mathrm{eff} (r)\) is the effective temperature in the Novikov–Thorne model defined as \(\mathcal {F}(r) = \sigma T^4_\mathrm{eff}\). \(\sigma \) is the Stefan–Boltzmann constant and \(\mathcal {F}(r)\) is the time-averaged energy flux from the surface of the disk,

$$\begin{aligned} \mathcal {F}(r) = \frac{\dot{M}}{4 \pi \sqrt{-G}} \frac{-\partial _r \Omega }{(E - \Omega L)^2} \int _{R_\mathrm{ISCO}}^{r} (E - \Omega L) (\partial _\rho L) \text {d}\rho ,\nonumber \\ \end{aligned}$$
(60)

\(\dot{M}\) is the mass accretion rate. E, L, and \(\Omega \) are, respectively, the specific energy, the axial component of the specific angular momentum, and the orbital frequency of equatorial circular orbits, while G is the determinant of the near equatorial plane metric.

Observations measure the flux at the distance \(r_\mathrm{o}\) of the observer. The impact of our deformation parameters on the thermal spectrum of thin disks is illustrated in Figs. 6 and 7. The results could have been imagined from the contour levels of the Novikov–Thorne radiative efficiency \(\eta _\mathrm{NT}\) [20]. The shape of the spectrum is simple, as it is just a multi-blackbody spectrum because the disk radiates as a blackbody locally and then one has to integrate radially. The high energy cut-off of the spectrum is determined by the Novikov–Thorne radiative efficiency, and therefore the measurement of the thermal component of the disk roughly corresponds to the measurement of \(\eta _\mathrm{NT}\). We cannot distinguish objects with the same \(\eta _\mathrm{NT}\) from the sole observation of the disk’s thermal spectrum, and therefore objects on the same contour level of \(\eta _\mathrm{NT}\) have substantially the same spectrum. We note that in our figures we have considered even objects with spin parameter larger than 1. In the Kerr metric, for \(|a_*| > 1\) there is no black hole and there are reasons to ignore these objects (see e.g. the discussion in Ref. [8]). In our case, these objects are black holes and they cannot be excluded a priori.

8 Concluding remarks

In this work, we have discussed a simple parametrization to describe possible deviations from the Kerr metric and test astrophysical black hole candidates. Our metric is suitable for the numerical calculations required to evaluate the electromagnetic spectrum of thin disks, in particular the thermal spectrum and the iron line profile. These are today the two leading techniques to probe the spacetime geometry around black hole candidates [24].

The continuum-fitting method is normally used for stellar-mass black hole candidates only, because the temperature of a Novikov–Thorne disk scales as \(M^{-0.25}\): for \(M \approx 10\) \(M_\odot \), the spectrum is in the soft X-ray band, while for supermassive black hole candidates it is in the UV/optical bands, where dust absorption makes an accurate measurement impossible. Stronger constraints on possible deviations from the Kerr metric can be obtained when the inner edge of the disk is very close to the compact object [43]. The best target may be the black hole binary Cygnus X-1 [44].

The iron line method can be used for both stellar-mass and supermassive black hole candidates, because the measurement does not directly depend on the mass of the object. Stellar-mass black hole candidates have the advantage to be brighter, so the photon count number in the iron line is high enough, which is not the case in AGN data. However, the spectrum of black hole binaries is more difficult to model (mainly because of the higher temperature of the disk) and the low energy tail of the iron line overlaps with the thermal component of the disk. At the moment it is not clear if the best targets to test the Kerr metric with the iron line are stellar-mass or supermassive black hole candidates. Even in this case, stronger constraints can be obtained if the inner part of the accretion disk is very close to the central object. The sources should also be sufficiently reflection-dominated, to have a stronger contrast between the primary and the reflection component. Cygnus X-1 should again be one of the most promising targets in the case of black hole binaries [45]. For AGN, current observations suggest that good targets would be NGC1365 [46] and 1H0707-495 [47], as both sources are bright, their reflection component is strong, and the inner edge of the disk seems to be very close to the black hole candidate.

In the analysis of real data, it is usually necessary to compute many spectra to fit the model parameters. In the Kerr case, it is possible to exploit a number of nice properties, with the results that numerical calculations can be fast and accurate. These properties are usually absent in non-Kerr metrics, and the calculation times become too long. To have a rough idea, the ray tracing calculations for one spectrum in the Kerr metric take about 5 min with a standard computer, so the calculation of a grid of, say, 40 spins and 20 angles takes something like 3 days. If the calculations are done by solving the geodesic equations, the computation time increases by about an order of magnitude is we want the same calculation accuracy. Moreover, we do not have a 2D grid any more, but a 3D grid if we consider one deformation parameter, or a 4D grid if we have two deformation parameters at the same time. Any parametrization has its own advantages and disadvantages, and in any case we cannot pretend to test black hole candidates with the most general black hole solution, because this would require an infinite number of deformation parameters and even in the presence of high quality data it is impossible to constrain several deformation parameters at the same time. Bearing in mind that any parametrization is necessary ad hoc and subject to criticisms, in our case we have looked for something suitable to analyze X-ray data.