# Note on equatorial geodesics in circular spacetimes

K. van Aelst

► To cite this version:

K. van Aelst. Note on equatorial geodesics in circular spacetimes. Classical and Quantum Gravity, 2020, 37 (20), pp.207001. 10.1088/1361-6382/aba80c. hal-02973190

HAL Id: hal-02973190

https://hal.science/hal-02973190v1

Submitted on 26 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Note on equatorial geodesics in circular spacetimes

K. Van Aelst¹

¹Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190 Meudon, France

E-mail: karim.van-aelst@obspm.fr

**Abstract.** General results on equatorial geodesics are exposed in the case of circular spacetimes featuring an equatorial reflection symmetry. The way the geodesic equation equivalently rewrites in terms of an effective potential is explicitly recalled for circular and non-circular equatorial geodesics. This provides a practical tool to discuss their stability. Statements are illustrated in Kerr spacetime.

**Keywords:** geodesics, stability, circular spacetime

## 1. Introduction

One of the first reference works on equatorial geodesics date back to the investigations of Boyer and Price [1], who focused on the Kerr metric just before the seminal articles of Carter on general Kerr geodesics [2, 3, 4]. Later on, Bardeen [5], Press and Teukolsky [6] considered equatorial geodesics in circular spacetimes. A defining property of the latter is to admit a quasi-isotropic coordinate system, i.e. a coordinate system $(t, r, \theta, \phi)$ with respect to which the line element writes

$$ds^2 = -N^2dt^2 + A^2(dr^2 + r^2d\theta^2) + B^2r^2\sin^2\theta(d\varphi - \omega dt)^2, \quad (1)$$

where the metric functions $N$, $A$, $B$ and $\omega$ only depend on coordinates $r$ and $\theta$.

At the time, these authors were interested in the precession of the periapsis of bounded orbits, stability of thin dust disks and other processes around rotating black holes. In regard of the highly accurate observations realized by instruments like GRAVITY [7, 8, 9] and the Event Horizon Telescope [10, 11, 12], such investigations are still essential today, e.g. to examine orbits of stars, accretion disks around compact objects and the images they produce. Contemporary discussions on these topics can be found in references [13, 14, 15, 16, 17].

Circular spacetimes (1) represent a large subclass of stationary and axisymmetric spacetimes (the corresponding Killing vectors being $\partial_t$ and $\partial_\phi$) possessing the additional symmetry $(t, \phi) \mapsto (-t, -\phi)$. The present note further assumes the equatorial reflection symmetry

$$\forall \mu, \nu, g_{\mu\nu}(r, \pi/2 - \theta) = g_{\mu\nu}(r, \theta). \quad (2)$$
Such metrics notably include the Kerr family‡, numerical black hole metrics [19, 20, 21], models of rotating neutron stars [22] and boson stars [15]. Yet, note for instance that Kerr-Newman-Taub-NUT metric [23, 24, 25] does not fall into this family as equatorial symmetry (2) cannot hold in presence of the gravitomagnetic parameter. One may consult [26, 27, 28, 13, 29] and references therein for further technical details on circular spacetimes and comments on their physical relevance.

## 2. Conservation equations

This section recalls the elements relevant in investigating the existence and properties of the trajectory of a free massive (resp. massless) particle in the equatorial plane of a circular spacetime. To effectively search for such a trajectory, its parametrization is set to be the only one whose corresponding tangent vector is the 4-momentum $p$ of the particle. More explicitly, one looks for a timelike (resp. null), future-oriented curve

$$ \mathcal{C}: \lambda \mapsto (x^{\mu}(\lambda)) = (t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)) \quad (3) $$

such that the 4-momentum of the particle is $p^{\mu} = \dot{x}^{\mu}$, where a dot denotes differentiation with respect to the parameter $\lambda$. In addition, the particle is free if and only if its 4-momentum is parallely transported along its trajectory $\mathcal{C}$:

$$ \nabla_{\mu} p^{\mu} = 0, \quad (4) $$

so that $\mathcal{C}$ is an affinely parametrized geodesic by definition.

The geodesic equation (4) implies that the mass

$$ m = \sqrt{-p^2} \quad \text{is conserved along } \mathcal{C}. \quad (5) $$

In particular, if $m > 0$, $\lambda$ is necessarily the curvilinear abscissa (i.e. proper time) $\tau$ along $\mathcal{C}$ divided by $m$.

In a stationary and axisymmetric spacetime such as (1), equation (4) also implies conservation of the Killing energy and angular momentum:

$$ E = -\partial_t \cdot p \quad \text{is conserved along } \mathcal{C}, \quad (6) $$

$$ L = \partial_{\phi} \cdot p \quad \text{is conserved along } \mathcal{C}. \quad (7) $$

Such quantities are actual observables only if the particle ever reaches spacelike infinity, where they are the energy and angular momentum effectively measured by a zero angular momentum observer (ZAMO)§. In this case, $E \ge 0$ necessarily.

Finally, the trajectory is requested to be equatorial:

$$ \theta = \pi/2 \quad \text{is conserved along } \mathcal{C}, \quad (8) $$

<footnote>‡ The corresponding expressions of the metric functions $N$, $A$, $B$ and $\omega$ (and the transformation from the usual Boyer-Lindquist coordinates to quasi-isotropic coordinates) are given in appendix Appendix A; reference [18] may also be consulted.</footnote>

<footnote>§ The ZAMO are characterized by a 4-velocity colinear to $\nabla t$; as a result, one may check that the ZAMO are not freely falling, yet they fulfill property (7) with $L = 0$, hence their name.</footnote>
which implies $p^{\theta} = 0$.

The conservation equations (5), (6), (7) and (8) are thus four necessary conditions for a curve $C$ to describe an equatorial trajectory of a free particle.

### 3. Non-circular geodesics

For non-circular orbits, i.e. for any trajectory such that $p^r \neq 0$ ¶ almost everywhere$^+$, these four conservation equations are sufficient: they imply the geodesic equation (4) for the following reasons. The Killing equation for $\partial_t$

$$[\nabla_{\mu} \partial_{t}]_{\nu} + [\nabla_{\nu} \partial_{t}]_{\mu} = 0 \qquad (9)$$

and the Killing energy conservation (6) establish the covariant $t$ component of the geodesic equation:

$$[\nabla_{p} p]_{t} = \partial_{t} \cdot \nabla_{p} p = \nabla_{p} (\partial_{t} \cdot p) = -\nabla_{p} E = 0. \qquad (10)$$

The analogous argument for $\partial_{\phi}$ and (7) yield $[\nabla_{p} p]_{\phi} = 0$.

Besides, the covariant $\theta$ component of the geodesic equation also vanishes:

$$[\nabla_{p} p]_{\theta} = A^{2} r^{2} [\nabla_{p} p]^{\theta} = A^{2} r^{2} p^{\mu} p^{\nu} \Gamma^{\theta}_{\mu\nu} = 0 \qquad (11)$$

since $\Gamma^{\theta}_{r\theta}$ gets multiplied by $p^{\theta} = 0$, and all the other Christoffel symbols $\Gamma^{\theta}_{\mu\nu}$ vanish in the equatorial plane as sums of terms that are proportional either to $\cos\theta$, or to some angular derivative of the metric $\partial_{\theta}g_{\mu\nu}$, which is necessarily zero under the natural assumption of equatorial symmetry (2). When the latter does not hold, e.g. in Kerr-Newman-Taub-NUT spacetime, the covariant $\theta$ component of the geodesic equation might not come out so simply: it may require alternative constraints on the metric components, or invoke some of the other conservation equations. It is indeed interesting to note that each of the above three covariant components of the geodesic equation do not require any of the two other conservation equations to be derived. As a result, any trajectory satisfying conservation of $E$ (resp. $L$, resp. $\theta = \pi/2$ when equatorial symmetry holds) always satisfies the covariant $t$ (resp. $\phi$, resp. $\theta$) component of the geodesic equation.

Finally, mass conservation (5) rewrites as

$$0 = \nabla_{p} (p \cdot p) = p^{\mu} [\nabla_{p} p]_{\mu} = p^{r} [\nabla_{p} p]_{r}, \qquad (12)$$

which implies $[\nabla_{p} p]_{r} = 0$ since $p^{r} \neq 0$ almost everywhere.

Before treating the circular case $p^r = 0$, it is very useful to note that equations (6), (7) and (8) allow to rewrite the mass conservation equation (5) as a familiar first order ordinary differential equation on the radial coordinate function $r$:

$$\frac{\dot{r}^2}{2} + \mathcal{V}(r, m, E, L) = 0, \qquad (13)$$

|| This is equivalent to $p_{\theta} \neq 0$ in quasi-isotropic coordinates (1).

¶ This is equivalent to $p_{r} \neq 0$ in quasi-isotropic coordinates (1).

+ Radial momentum $p^r$ may only cancel at periapsis and apoapsis, when they exist.
where the effective potential $\mathcal{V}$ is defined as

$$ \mathcal{V}(r, m, E, L) = \frac{1}{2A^2} \left[ m^2 - \left( \frac{E - \omega L}{N} \right)^2 + \left( \frac{L}{Br} \right)^2 \right]. \qquad (14) $$

Since $\lambda = \tau/m$ for massive particles, note that (13) rewrites as

$$ \frac{1}{2} \left( \frac{dr}{d\tau} \right)^2 + \mathcal{V}(r, 1, \bar{E}, \bar{L}) = 0, \qquad (15) $$

where $\bar{E} = E/m$ and $\bar{L} = L/m$, so that the trajectories of free massive particles only depend on their Killing energy and angular momentum per unit mass.

Based on the explicit form (13) of the mass conservation equation, conditions on $\mathcal{V}$ and its partial derivatives will also allow to characterize the circular geodesics (see section 4) and study their stability (see section 5). So far, simply note that non-circular geodesics necessarily satisfy $\mathcal{V} < 0$ almost everywhere, while circular geodesics necessarily satisfy $\mathcal{V} = 0$ everywhere.

Finally, the procedure to explicitly construct all non-circular equatorial trajectories of free particles is to first pick an initial radial coordinate $r_0$, a Killing energy $E$, a Killing angular momentum $L$ and a mass $m$ such that $\mathcal{V}(r_0, m, E, L) < 0$ and $E - \omega(r_0)L > 0$ to guarantee that the trajectory is initially causal future-oriented (see equation (17) below); in particular, this necessarily requires $E > 0$ if $r_0$ is outside the ergoregion. Recall that the latter is the domain over which the pseudo-stationary Killing vector $\partial_t$ is non-timelike, i.e. $N^2 \le (\omega Br \sin\theta)^2$. Thus, realistic observers can no longer have constant spatial coordinates. Theoretically, this region allows to use particles to extract rotational energy from a black hole, as described by the Penrose process. In practice, the latter is not efficient enough to be significantly involved in astrophysical processes such as the relativistic jets emerging e.g. from quasars, although this used to be conjectured. Yet higher efficiencies might be reached around other objects, such as naked singularities or wormholes, or through more elaborate avatars of the process, such as the collisional Penrose process or superradiance. One may consult [30, 31] and references therein for discussions of these topics.

Secondly, the right-hand side of

$$ \dot{r} = \pm \sqrt{-2\mathcal{V}(r, m, E, L)} \qquad (16) $$

is sufficiently regular for equation (16) to admit a unique solution $\lambda \mapsto r_s(\lambda)$ once the $\pm$ sign is chosen (to determine whether the initial direction is ingoing or outgoing). Conservation of $E$ and $L$ then provide the solutions for $t$ and $\phi$:

$$ t_s(\lambda) = \int \frac{E - \omega(r_s)L}{N(r_s)^2} d\lambda, \qquad (17) $$

$$ \phi_s(\lambda) = \int \left[ \frac{L}{(B(r_s)r_s)^2} + \omega(r_s) \frac{E - \omega(r_s)L}{N(r_s)^2} \right] d\lambda, \qquad (18) $$

where the second argument $\theta_s = \pi/2$ of the metric functions is omitted without loss of clarity. Note that frame dragging can be read off from relation (18): radial free fall from infinity (hence $L = 0$) does not remain radial in a rotating spacetime as $\phi_s$ receives contribution from the non-vanishing metric function $\omega$.
## 4. Circular geodesics

As mentioned earlier, $p^r = 0$ (circular orbit) forbids to establish the radial geodesic equation from mass conservation, which is then redundant as a linear combination of the $E$ and $L$ conservation equations (i.e. the covariant $t$ and $\phi$ geodesic equations). Therefore, one additional equation is missing to realize a geodesic. Indeed, there are so far multiple solutions to the problem {$(6), (7), (8), p^r = 0$} (which implies mass conservation): for any $E, L$ and $r_0$, the curve

$$ \lambda \mapsto \left( \frac{E - \omega_0 L}{N_0^2} \lambda, r_0, \pi/2, \left[ \frac{L}{(B_0 r_0)^2} + \omega_0 \frac{E - \omega_0 L}{N_0^2} \right] \lambda \right) \quad (19) $$

is circular (an index 0 means that the metric function is evaluated at $r_0$, in the equatorial plane) with conserved Killing energy $E$ and angular momentum $L$. Additionally, one may simply require $E \ge \omega_0 L + N_0|L|/(B_0 r_0)$ (which is always positive outside the ergoregion) to describe a causal future-oriented trajectory; this amounts to requiring that mass conservation $V(r_0, m, E, L) = 0$ holds for a real constant $m$ (i.e. $m^2 > 0$).

All these circular orbits are distinct solutions to the same problem, but at most one of them is a geodesic, sometimes none (intuitively, the other orbits are accelerated inward if they rotate “faster” than a geodesic to compensate for the centrifugal effect, outward otherwise). To search for a geodesic among them, one obtains an additional prescription from differentiating equation (13):

$$ \dot{r} [\ddot{r} + V'(r, m, E, L)] = 0, \quad (20) $$

where $'$ denotes differentiation with respect to the first argument (the radial coordinate).

Of course, for any circular orbit (19), equation (20) holds because $\dot{r} = 0$, while it implies

$$ \ddot{r} + V'(r, m, E, L) = 0, \quad (21) $$

for all non-circular geodesics. The missing condition to realize a circular geodesic at $r_0$ is then obtained by requiring equation (21) to hold even in the circular limit, i.e. when $\ddot{r} = 0$, yielding

$$ V'(r_0, m, E, L) = 0. \quad (22) $$

On figure 1a for instance, circular geodesics exist at the zeros of the bottom red and top blue curves since the latter cancel in a stationary way (the corresponding geodesics are respectively represented by the red and blue dashed lines in figure 1b), while the zeros of the two other curves can only correspond to accelerated circular orbits or to the periapsis and apoapsis of a non-circular geodesic; more precisely, the smallest (resp. greatest) zero of the top red curve is the unique apoapsis (resp. periapsis) of a geodesic reaching the event horizon (resp. infinity, as represented by the red solid line in figure 1b) while the smallest (resp. greatest) zero of the bottom blue curve is the periapsis (resp. apoapsis) of a bounded geodesic (represented by the blue solid line in figure 1b). Note here that, in the asymptotically flat case, $N \to 1$, $A \to 1$, $B \to 1$ and $\omega \to 0$ at infinity, so that $V(r, m, E, L) \to (1 - E^2)/2$; therefore, a non-circular unbounded geodesic (such
as the one represented by the red solid line on figure 1b ruled by the right negative branch of the top red curve of figure 1a) requires $E > 1$.

To show that the additional condition (22) allows to establish the $r$ geodesic equation, define

$$W(r, m, E, L) = m^2 - \left(\frac{E - \omega L}{N}\right)^2 + \left(\frac{L}{Br}\right)^2, \qquad (23)$$

so that $V(r, m, E, L) = W(r, m, E, L)/(2A^2)$, and hence

$$V'(r_0, m, E, L) = \left(\frac{1}{2A^2}\right)'_0 W(r_0, m, E, L) + \frac{W'(r_0, m, E, L)}{2A_0^2} \qquad (24)$$

$$= \frac{W'(r_0, m, E, L)}{2A_0^2}, \qquad (25)$$

since mass conservation $V(r_0, m, E, L) = 0$ is equivalent to $W(r_0, m, E, L) = 0$. Therefore, condition (22) is equivalent to

$$0 = \frac{W'(r_0, m, E, L)}{2} \qquad (26)$$

$$= \frac{E - \omega_0 L}{N_0^2} \left[ (E - \omega_0 L) \frac{N_0'}{N_0} + L \omega_0' \right] - \left( \frac{L}{B_0 r_0} \right)^2 \left( \frac{B_0'}{B_0} + \frac{1}{r_0} \right). \qquad (27)$$

Then, injecting the explicit expressions

$$L = B_0^2 r_0^2 (p^\phi - \omega_0 p^t), \qquad (28)$$

$$E = N_0^2 p^t + \omega_0 L, \qquad (29)$$

into (26) immediately yields

$$0 = p^t [N_0 N_0' p^t + B_0^2 r_0^2 (p^\phi - \omega_0 p^t) \omega_0'] - B_0^2 r_0^2 (p^\phi - \omega_0 p^t)^2 \left(\frac{B_0'}{B_0} + \frac{1}{r_0}\right) \qquad (30)$$

$$= A_0^2 \left[ \Gamma^r_{tt} (p^t)^2 + 2 \Gamma^r_{t\phi} p^t p^\phi + \Gamma^r_{\phi\phi} (p^\phi)^2 \right] \qquad (31)$$

$$= A_0^2 \Gamma^r_{\mu\nu} p^\mu p^\nu = A_0^2 [\nabla_p p]^r = [\nabla_p p]_r. \qquad (32)$$

Therefore, circular geodesics are precisely the circular orbits (19) that continue property (21) of non-circular geodesics. Gathered with the result of section 3, the circular and non-circular free trajectories are characterized as follows. A curve $C: \lambda \mapsto (x^\mu(\lambda))$ describes an equatorial free particle with 4-momentum $p^\mu = \dot{x}^\mu$ if and only if it satisfies the three conservation equations $\{(6), (7), (8)\}$ and either one of the two following conditions:

* the fourth conservation equation (5) and $V \neq 0$ almost everywhere on $C$ (in this case, $C$ is non-circular);

* $V = 0$ and $V' = 0$ everywhere on $C$ (in this case, $C$ is circular).

Finally, to explicitly construct all free circular equatorial trajectories, recall that all the curves (19) satisfying $E \ge \omega L + N|L|/(Br)$ (the indices 0 are now removed although all the statements in the remaining of the section will only apply to circular orbits) are very good candidates because they satisfy the three conservation equations $\{(6), (7), (8)\}$
<figure><img src="image_2.png" /><figcaption>(a) Effective potentials of timelike orbits with respect to $\bar{r} = r/r_H$ for different Killing energy and angular momentum per unit mass (only the ratios $\bar{E} = E/m$ and $\bar{L} = L/m$ are relevant in the massive case). The bottom red (resp. top blue) curve corresponds to the Killing energy and angular momentum of the timelike circular geodesic at $4r_H$ (resp. $15r_H$) while the other red (resp. blue) curve has same $\bar{L}$ but a slightly smaller (resp. greater) $\bar{E}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) The red (resp. blue) dashed circle corresponds to the zero of the bottom red (resp. top blue) curve of figure 1a. The red (resp. blue) solid line is ruled by the right, negative, unbounded (resp. bounded) branch of the top red (resp. bottom blue) curve of figure 1a. Plot realized with the ray-tracing code *GYOTO* [32].</figcaption></figure>

<caption>Figure 1: Potentials and spatial projections of timelike circular and non-circular equatorial geodesics in a Kerr spacetime ($a/M \approx 0.52$).</caption>
together with $V(r, m, E, L) = 0$, since this relation is the definition of $m$ for a circular orbit at $r$. It thus only remains to derive which final constraint emerges from requiring $V'(r, m, E, L) = 0$ (in which $m$ actually does not appear because of $V(r, m, E, L) = 0$). This constraint will first be formulated in terms of the “signed norm” $V$ of the spatial velocity $v$ measured by the ZAMO, defined below. From $V$ will then be deduced the values of $E$ and $L$ to be injected into (19) to define a free circular trajectory at $r$.

For a particle with mass $m$, first denote $\mathcal{E}$ and $v$ the energy and spatial velocity measured by the ZAMO*, and $n$ the 4-velocity of the latter, so that the 4-momentum of the particle decomposes as

$$p = \mathcal{E}(n + v) \text{ with } n \cdot v = 0. \qquad (33)$$

One obtains

$$\mathcal{E} = \frac{E - \omega L}{N} \qquad (34)$$

and

$$v = \frac{V}{Br} \partial_{\phi} \qquad (35)$$

where

$$V = \frac{L}{Br\mathcal{E}} \qquad (36)$$

Note that $v^2 = V^2$, hence the name “signed norm” for $V$. Also recall here that a trajectory is defined to be prograde (resp. retrograde) when $\omega L > 0$ (resp. $\omega L < 0$); since $\omega$ does not generally cancel, the $\phi$ coordinate may be chosen such that $\omega > 0$, which then simplifies in the definition. Based on relations (35) and (36), and the fact that $\omega(r_H)$ equals the angular velocity $\Omega_H$ of the event horizon‡, a prograde (resp. retrograde) trajectory intuitively rotates in the same (resp. opposite) direction as the black hole for the ZAMO. However, for an observer at infinity, the angular velocity $p^\phi/p^t$ of a circular orbit (19) is

$$\frac{p^{\phi}}{p^{t}} = \omega + \frac{NL}{Br\sqrt{(mBr)^{2} + L^{2}}}, \qquad (37)$$

so that one of the retrograde orbit may appear prograde from infinity if and only if it is outside the ergoregion and $L > -m\omega(Br)^2/\sqrt{-N^2 + (\omega Br)^2}$.

Note also that equation (26) (which is equivalent to $V' = 0$ as $V = 0$) is independent of $m$ and hence homogeneous with respect to $E$ and $L$. Therefore, injecting relations (34) and (36) into this expression allows to simplify all $\mathcal{E}$, which finally yields the following second order equation in $V$:

$$\left(\frac{B'}{B} + \frac{1}{r}\right) V^2 - \frac{Br\omega'}{N} V - \frac{N'}{N} = 0. \qquad (38)$$

<footnote>* For a massive particle, $\mathcal{E} = \Gamma m$ where $\Gamma$ is the Lorentz factor of the particle with respect to the ZAMO; for a massless particle, $\mathcal{E} = h\nu$ where $\nu$ is the frequency measured by the ZAMO.</footnote>

<footnote>‡ In quasi-isotropic coordinates, the event horizon is always located at a constant radial coordinate $r_H$.</footnote>
As announced, only the spatial velocity of the particle is constrained and not $\mathcal{E}$, which means that, where a circular timelike geodesic exists, it can be the worldline of any massive particle regardless of its mass provided it has the right velocity, and where a circular null geodesic exists (i.e. a photon ring), it can be the worldline of any photon regardless of its frequency (this all seems consistent with the equivalence principle). For a timelike circular geodesic to exist at $r$, the values at $r$ of the metric functions involved in (38) need to be such that at least one of the roots, if any exists, belongs to $(-1, 1)$ (the ZAMO must measure subluminal velocities). In the massless case, $V = \pm 1$ so that at any photon ring, if any exists, the metric functions need to be such that 1 or $-1$ is a root of equation (38). One then only needs to study the roots $V_{\pm}$ of equation (38) to conclude about existence and location of timelike circular geodesics and photon rings. These roots exist if and only if the discriminant

$$D = \left(\frac{Br\omega'}{N}\right)^2 + \frac{4N'}{N}\left(\frac{B'}{B} + \frac{1}{r}\right) \qquad (39)$$

of equation (38) is non-negative, in which case one has

$$V_{\pm}(r) = \frac{\frac{Br\omega'}{N} \pm \sqrt{D}}{2\left(\frac{B'}{B} + \frac{1}{r}\right)}. \qquad (40)$$

In Kerr spacetime, for any angular velocity $\Omega_H^{\dagger\dagger}$, function $D$ monotonically decreases from an infinite value at the horizon down to a zero limit as $r \to +\infty$ (see figure 2a). Therefore, each velocity function $V_{\pm}$ necessarily becomes luminal at some point corresponding to a photon ring (marked with a vertical line from 0 to 1 in figure 2b), beyond which timelike circular geodesics exist everywhere.

Finally, relations (34) and (36) then yield

$$L_{\pm} = \mathcal{E}BrV_{\pm}, \qquad (41)$$

$$E_{\pm} = \mathcal{E}(N + Br\omega V_{\pm}), \qquad (42)$$

i.e. in the massive case,

$$L_{\pm} = \Gamma_{\pm}mBrV_{\pm}, \qquad (43)$$

$$E_{\pm} = \Gamma_{\pm}m(N + Br\omega V_{\pm}), \qquad (44)$$

with

$$\Gamma_{\pm} = (1 - V_{\pm}^2)^{-1/2}, \qquad (45)$$

while in the massless case,

$$L_{\pm} = \pm h\nu Br, \qquad (46)$$

$$E_{\pm} = h\nu(N \pm Br\omega). \qquad (47)$$

These quantities are plotted on figure 2c, 2d and 2e in the massive case for different dimensionless angular velocities $\bar{\Omega}_H = r_H\Omega_H$ in Kerr spacetime. In particular, the “+”

$^{\dagger\dagger}$In Kerr spacetime, angular velocity $\Omega_H$ relates the usual mass $M$ and spin parameter $a$ according to $\Omega_H = a/(2M(M + \sqrt{M^2 - a^2}))$.
<figure><img src="image_4.png" /><figcaption>(a) Positivity of *D* allows (possibly superluminal) circular geodesics.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Velocities and photon rings.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Lorentz factor of timelike circular geodesics.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Killing energy of timelike circular geodesics.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) Killing angular momentum of timelike circular geodesics.</figcaption></figure>

<caption>Figure 2: Kinematic characteristics of circular geodesics in Kerr spacetime for different dimensionless angular velocities $\bar{\Omega}_{\mathcal{H}} = r_{\mathcal{H}}\Omega_{\mathcal{H}}$.</caption>
and “-” quantities are no longer merely equal or opposite in the rotating cases and respectively follow the evolution of the “+” and “-” photon rings in figure 2b: the “+” (resp. “-”) Lorentz factor and Killing energy are for instance roughly shifted to the left (resp. right) of the Schwarzschild profile.

## 5. Stability of circular geodesics

For circular geodesics, the radial equation (13) expectedly provides a stability criteria based on convexity. A non-constant perturbation $\delta$ (a constant perturbation would not threaten stability) to a circular geodesic at $r$ allows to use equation (21) instead:

$$ \ddot{\delta} + \mathcal{V}'(r + \delta, m, E + \delta_E, L + \delta_L) = 0 \qquad (48) $$

in which Taylor expanding and invoking condition (22) yields

$$ \ddot{\delta} + \mathcal{V}''(r, m, E, L)\delta = O(\delta_E) + O(\delta_L) + O(\delta^2 + \delta_E^2 + \delta_L^2). \qquad (49) $$

If $\mathcal{V}''(r, m, E, L) < 0$, then $\delta$ must be accelerating away from $r$ to preserve the asymptotic orders in the right-hand side. To guarantee that any $\delta$ is bounded in some neighbourhood of $r$ then requires positive $\mathcal{V}''(r, m, E, L)$.

Actually, the values of $E$ and $L$ for a circular geodesic at $r$ are necessarily $E_{\pm}(r)$ and $L_{\pm}(r)$, explicitly given by relations (41) and (42). In practice, one should thus study the sign of the two functions

$$ \mathcal{V}_{\pm}'' : r \mapsto \mathcal{V}''(r, m, E_{\pm}(r), L_{\pm}(r)), \qquad (50) $$

on the set on which the discriminant $D$ is non-negative. Actually, the expressions $\mathcal{V}''(r, m, E_{\pm}(r), L_{\pm}(r))$ are homogeneous with respect to $\mathcal{E}$, so that their sign do not depend on $\Gamma_{\pm}m$ in the massive case nor on $h\nu$ in the massless case. Therefore, the stability of causal circular geodesics only depends on the sign of the two functions (50) and concerns massive particles where $V_{\pm}(r) \in (-1, 1)$ and massless ones where $V_{\pm}(r) = \pm 1$, regardless of whether the expressions used for $E_{\pm}$ and $L_{\pm}$ apply to a massive or a massless particle.

Figures 3a and 3b illustrate all this in the Schwarzschild case: for each radii $r_0$ marked with a dashed vertical line of a given colour in figure 3a, the potential $\mathcal{V}(\cdot, m, E_{\pm}(r_0), L_{\pm}(r_0))$, which corresponds to the geodesic at $r_0$, is plotted with the same colour, and thus cancels in a stationary way at $r_0$. Then, for each $r_0$, function $\mathcal{V}_{\pm}''$ in figure 3b extracts the convexity of $\mathcal{V}(\cdot, m, E_{\pm}(r_0), L_{\pm}(r_0))$ at the corresponding $r_0$: the purple curve is concave at $r_0 = 5r_{\mathcal{H}}$ on figure 3a (hence unstable circular geodesic), so that $\mathcal{V}_{\pm}''$ is negative at $5r_{\mathcal{H}}$ on figure 3b, whereas the turquoise curve is convex at $r_0 = 20r_{\mathcal{H}}$, so that $\mathcal{V}_{\pm}''$ is positive at $20r_{\mathcal{H}}$. The limiting case (grey lines) such that $\mathcal{V}(\cdot, m, E_{\pm}(r_0), L_{\pm}(r_0))$ cancels as an inflection point at $r_0 \simeq 9.9r_{\mathcal{H}}$ ($\mathcal{V}_{\pm}''(r_0) = 0$) defines the innermost stable circular orbit (ISCO)$^\dagger$.

<footnote>$^\dagger$ In Boyer-Lindquist coordinates, the ISCO radius $R_0$ of Schwarzschild spacetime is known to equal $6M$; injecting this value into equation (A.8) does yield $r_0 = (5 + 2\sqrt{6})M/2 \simeq 9.9r_{\mathcal{H}}$.</footnote>
<figure><img src="image_9.png" /><figcaption>(a) Potentials $V(\cdot, m, E_\pm(r_0), L_\pm(r_0))$ of circular(b) $V''_\pm(r_0)$ is the convexity of $V(\cdot, m, E_\pm(r_0), L_\pm(r_0))$ geodesics at various $r_0$ in Schwarzschild spacetime. evaluated at $r_0$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 3: Positivity of $V''_\pm$ rules stability and hence location of ISCO.</figcaption></figure>

Finally, figure 3c gathers the functions $V''_\pm$ for different angular velocities. The functions $V'_+$ and $V''_-$ are no longer equal in the rotating case ($V''_+$ globally increases with rotation while $V''_-$ globally decreases) and thus respectively define an ISCO. Based on figure 2b, both ISCO are always located beyond the corresponding photon ring, so that the latter are always unstable.

Note that for a stable circular geodesic ($r_0, m, E, L$), the function $V(\cdot, m, E, L)$ realizes a local minimum (equal to zero) at $r_0$, so that it is strictly positive in a neighbourhood of $r_0$ except at $r_0$. Yet, one always has

$$ \partial_E V(r_0, m, E, L) = -\frac{\mathcal{E}}{N_0 A_0^2} < 0 \qquad (51) $$

since $\mathcal{E}$ is always strictly positive for a causal future-oriented curve. Decreasing $E$ thus
<figure><img src="image_11.png" /><figcaption>Figure 4: Effective potential $V^{(\bar{L})}$ (yellow surface) of massive particles at fixed $\bar{L}$ in a Kerr spacetime ($a/M \simeq 0.52$). Its zeros lie on the black curve $E_{\text{min}}^{(\bar{L})}$. Its profile at fixed Killing energy per unit mass $e$ (resp. $e'$) is highlighted as a red (resp. blue) curve on the blue (resp. red) vertical plane. $\bar{L}$ being fixed, the only circular geodesics are marked with the red and blue dots, where both the two curves to which they respectively belong are stationary, and whose common convexity determines their stability.</figcaption></figure>

increases $V$ locally‡, so that $V$ becomes strictly positive on a neighbourhood of $r_0$. But $V$ is necessarily negative or zero on any orbit ruled by (13), which means that, $m$ and $L$ being fixed, there can be no geodesic close to $(r_0, m, E, L)$ with smaller $E$: a stable circular geodesic at $r_0$ realizes a local minimum of $E$ on the set of geodesics having same $m$ and $L$.

As an additional note, fixing $m$ and $L$ also provides other interesting criteria to characterize circular geodesics among circular orbits (instead of $V' = 0$), and investigate their stability (instead of $V'' > 0$). For illustrative purposes, $V^{(m,L)}: (r, E) \mapsto V(r, m, E, L)$ is plotted for a Kerr spacetime in figure 4.

From this point of view, circular geodesics lie on the set of points $(r, E)$ such that $V^{(m,L)}(r, E) = 0$ (black horizontal curve in figure 4): it is the intersection of the image of $V^{(m,L)}$ (yellow surface) with the $V=0$ plane (green horizontal plane). For any of these points $(r_0, E_0)$ to actually correspond to a circular geodesic (rather than

<footnote>‡ This is illustrated on figure 1a: $E$ is decreased to switch from the bottom curve to the top curve of a given color.</footnote>
an accelerated circular orbit), condition (22) must hold, meaning that $r_0$ must be a stationary point of $V^{(m,L)}(\cdot, E_0)$, whose graph is the intersection of the image of $V^{(m,L)}$ (yellow surface) with the $E = E_0$ plane (blue or red vertical plane).

However, the set of points $(r, E)$ such that $V^{(m,L)}(r, E) = 0$ (black horizontal curve) is always the graph of a function well-defined with respect to $r$: $r \mapsto E(r)$. This comes from the fact that $V^{(m,L)}$ is a second order polynomial with respect to $E$. For any $r$, the function $V^{(m,L)}(r, \cdot)$ always admits two distinct roots:

$$E_{\text{min}}^{(m,L)}(r) = \omega(r)L + N(r)\sqrt{m^2 + \left(\frac{L}{B(r)r}\right)^2}, \quad (52)$$

$$E_{\text{neg}}^{(m,L)}(r) = -E_{\text{min}}^{(m,-L)}(r). \quad (53)$$

Since $V^{(m,L)}(r, \cdot)$ is a second order polynomial with negative dominant coefficient with respect to $E$, and since $E_{\text{min}}^{(m,L)}$ is always the greatest root, one has

$$\partial_E V^{(m,L)} (r, E_{\text{min}}^{(m,L)}(r)) < 0, \quad (54)$$

$$\partial_E V^{(m,L)} (r, E_{\text{neg}}^{(m,L)}(r)) > 0. \quad (55)$$

Based on equation (51), $E_{\text{neg}}^{(m,L)}$ can in no case correspond to a causal future-oriented curve. Actually, it merely corresponds to all the causal past-oriented circular and non-circular equatorial geodesics. This had to be expected since the effective potential $V$ and the three other conservation equations can as well be used to describe them. More precisely, based on relations (14), (16), (17) and (18), the past-oriented version of a future-oriented non-circular geodesic is obtained by switching the sign of both $E$ and $L$ and taking the opposite sign in (16). This is even easier to check in the circular case (19), and relation (53) confirms that $E_{\text{neg}}^{(m,L)}$ is merely the Killing energy of the past-oriented version of the circular orbit having opposite Killing angular momentum. One may thus focus on $E_{\text{min}}^{(m,L)}$ alone, which can only be negative for retrograde orbits inside the ergoregion.

By definition, one has

$$V^{(m,L)} (r, E_{\text{min}}^{(m,L)}(r)) = 0, \quad (56)$$

so that

$$E_{\text{min}}^{(m,L)'}(r) = -\frac{V^{(m,L)'}(r, E_{\text{min}}^{(m,L)}(r))}{\partial_E V^{(m,L)}(r, E_{\text{min}}^{(m,L)}(r))} \quad (57)$$

where division is allowed by (54). Therefore, the circular orbit $(r, m, E_{\text{min}}^{(m,L)}(r), L)$ is a geodesic if and only if $E_{\text{min}}^{(m,L)'}(r) = 0$ (instead of $V'(r, m, E_{\text{min}}^{(m,L)}(r), L) = 0$). Graphically, one may check on figure 4 that at the circular geodesic marked with a red (resp. blue) dot, both the red curve $V^{(\bar{L})}(\cdot, e)$ (resp. blue curve $V^{(\bar{L})}(\cdot, e')$) and function $E_{\text{min}}^{(\bar{L})}$ are stationary.
<figure><img src="image_12.png" /><figcaption>Figure 5: Function $\bar{E}_{min}^{(\bar{L})}$ of massive particles for various positive $\bar{L}$ in a Kerr spacetime ($a/M \simeq 0.52$).</figcaption></figure>

The second derivative of relation (56) at a circular geodesic ($r_0, m, E_{min}^{(m,L)}(r_0), L$) (to use $E_{min}^{(m,L)'}(r_0) = 0$) merely yields

$$E_{min}^{(m,L)''}(r_0) = -\frac{\mathcal{V}^{(m,L)''}(r_0, E_{min}^{(m,L)}(r_0))}{\partial_E \mathcal{V}^{(m,L)}(r_0, E_{min}^{(m,L)}(r_0))} \quad (58)$$

which has the same sign as $\mathcal{V}^{(m,L)''}(r_0, E_{min}^{(m,L)}(r_0))$ because of (54). Therefore, such circular geodesic is stable if and only if $E_{min}^{(m,L)''}(r_0) > 0$ (instead of $\mathcal{V}''(r_0, m, E_{min}^{(m,L)}(r_0), L) > 0$). Graphically, one may check on figure 4 that at the stable (resp. unstable) circular geodesic marked with a red (resp. blue) dot, both the red curve $\mathcal{V}^{(\bar{L})}(\cdot, e)$ (resp. blue curve $\mathcal{V}^{(\bar{L})}(\cdot, e')$) and function $E_{min}^{(\bar{L})}$ are convex (resp. concave).

It is thus possible to determine the location and stability of circular geodesics at fixed $m$ and $L$ from the function $E_{min}^{(m,L)}$ rather than $\mathcal{V}^{(m,L)}$. In this case, one would then study the dependence of $E_{min}^{(m,L)}$ on the parameter $L$ (graphically, changing $L$ deforms the yellow surface in figure 4 and hence the black curve). Figure 5 illustrates this for a Kerr spacetime. At fixed $L = L_0$, circular geodesics correspond to the stationary points. But then, for any stationary point ($r_0, E_0$), $E_0$ is necessarily equal to $E_{\pm}(r_0)$ (the sign depending on whether the orbit is prograde or retrograde) which is why the black dashed curve intersects each coloured curve at its stationary points.

Non-circular geodesics with Killing energy $E$ are possible where $E > E_{min}^{(m,L)}$.
Graphically, such a geodesic covers the region where the horizontal line E is above $E_{min}^{(m,L)}$ while the abscissae of their intersections locate the periapsis and apoapsis of the geodesic.

Despite providing another interesting point of view on geodesics and stability criteria, using $E_{min}^{(m,L)}$ to investigate stability is laborious since it requires to locate the circular geodesics and evaluate convexity for each L. Instead, it is much more efficient and exhaustive to focus on the sign of the two functions $V''_{\pm}$ given by (50), as illustrated for Kerr spacetime with figure 3c above.

For practical use, it is finally interesting to mention that the existence of a bounded non-circular orbit strongly suggests the existence of a stable circular orbit at some radius between its apsides, as one would expect intuitively. Denote $m, E$ and $L$ the mass, Killing energy and angular momentum of a particle following a bounded non-circular orbit such as the solid blue line in figure 1b ruled by the right, negative, bounded branch of the bottom blue curve of figure 1a. The argument is that $V(\cdot, m, E, L)$ necessarily admits a local minimum at some $r_0$ between its apsides, i.e. such that

$$V(r_0, m, E, L) < 0, \qquad (59)$$

$$V'(r_0, m, E, L) = 0, \qquad (60)$$

$$V''(r_0, m, E, L) > 0. \qquad (61)$$

Since $V(r, m, 0, 0) > 0$ and $V(r, m, E, L) < 0$, there necessarily exists a strictly positive factor $\alpha$ rescaling $E$ and $L$ in such a way that

$$V(r, m, \alpha E, \alpha L) = 0. \qquad (62)$$

Defining

$$\chi(r, m, E, L) = \frac{1}{2A^2} \left[ - \left( \frac{E - \omega L}{N} \right)^2 + \left( \frac{L}{Br} \right)^2 \right] \qquad (63)$$

so that $V = 1/(2A^2) + \chi$, one has, for any $\beta$,

$$V'(r_0, m, \beta E, \beta L) \approx \chi'(r_0, m, \beta E, \beta L) = \beta^2 \chi'(r_0, m, E, L), \qquad (64)$$

$$V''(r_0, m, \beta E, \beta L) \approx \chi''(r_0, m, \beta E, \beta L) = \beta^2 \chi''(r_0, m, E, L) \qquad (65)$$

far from the strong-field region (so that $A$ slowly varies) where most stars orbiting a central object would be observed. The quantities (64) and (65) cancel for $\beta = 1$ according to (60) and (61), which implies the same for $\beta = \alpha$. The parameters $(r_0, m, \alpha E, \alpha L)$ are thus close to define a stable circular orbit, and it should generally be possible to finish tuning them to obtain an exact stable circular orbit. This makes the previously discussed stability criteria for circular orbits relevant for observations although no exactly circular orbit exists; this means for instance that observing a star on a non-circular orbit around a black hole logically requires the ISCO to be located below the apoapsis of the star orbit.

## Acknowledgements

The author thanks Eric Gourgoulhon for helpful discussions and comments on this paper. This work was supported by the CNRS project 80PRIME-TNENGRAV.
Appendix A. Kerr metric in quasi-isotropic coordinates

Denoting *M* the mass and *a* the spin parameter, the four metric functions involved in the quasi-isotropic expression (1) of the Kerr metric explicitly write as

$$N^2 = \frac{\Sigma\Delta}{\Sigma(R^2 + a^2) + 2a^2MR\sin^2\theta}, \qquad (A.1)$$

$$A^2 = \frac{\Sigma}{r^2}, \qquad (A.2)$$

$$B^2 = \frac{1}{r^2} \left( R^2 + a^2 + \frac{2a^2MR\sin^2\theta}{\Sigma} \right), \qquad (A.3)$$

$$\omega = \frac{2aMR}{\Sigma(R^2 + a^2) + 2a^2MR\sin^2\theta}, \qquad (A.4)$$

where

$$R = r + \frac{M^2 - a^2}{4r} + M, \qquad (A.5)$$

$$\Sigma = R^2 + a^2 \cos^2 \theta, \qquad (A.6)$$

$$\Delta = R^2 + a^2 - 2MR. \qquad (A.7)$$

The function *R* defined by relation (A.5) actually corresponds to the radial coordinate of the Boyer-Lindquist system (which is presented e.g. in section 33.2 of [33]). It is inverted as

$$r = \frac{1}{2} \left( R + \sqrt{R^2 - 2MR + a^2} - M \right). \qquad (A.8)$$

The remaining coordinates are identical in the quasi-isotropic and Boyer-Lindquist systems.

References

[1] Boyer R H and Price T G 1965 *Math. Proc. Camb. Philos. Soc.* **61** 531-534

[2] Carter B 1966 *Phys. Rev.* **141** 1242-1247 URL https://link.aps.org/doi/10.1103/PhysRev.141.1242

[3] Carter B 1968 *Phys. Rev.* **174** 1559-1571 URL https://link.aps.org/doi/10.1103/PhysRev.174.1559

[4] Carter B 1968 *Commun. Math. Phys.* **10** 280-310 URL https://projecteuclid.org:443/euclid.cmp/1103841118

[5] Bardeen J 1970 *ApJ* **161** 103-9

[6] Bardeen J, Press W and Teukolsky S 1972 *ApJ* **178** 347

[7] Abuter R *et al.* 2017 *A&A* **602** A94 URL https://doi.org/10.1051/0004-6361/201730838

[8] GRAVITY Collaboration 2018 *A&A* **615** L15 URL https://doi.org/10.1051/0004-6361/201833718

[9] GRAVITY Collaboration 2018 *A&A* **618** L10 URL https://doi.org/10.1051/0004-6361/201834294

[10] Doleman S *et al.* 2009 *ASTRO2010 Decadal Review Panels* (Preprint 0906.3899)

[11] Issaoun S *et al.* 2019 *ApJ* **871** 30 URL https://iopscience.iop.org/article/10.3847/1538-4357/aaf732

[12] Collaboration T E *et al.* 2019 *ApJL* **875** 1 URL https://iopscience.iop.org/article/10.3847/2041-8213/ab0ec7
[13] Gourgoulhon E 2010 An introduction to the theory of rotating relativistic stars *CompStar 2010: School and Workshop on Computational Tools for Compact Star Astrophysics* Ganil, Caen, France, February 8-16, 2010 (Preprint 1003.5015)

[14] Pugliese D, Quevedo H and Ruffini R 2011 *Phys. Rev. D* **84** 044030 URL https://link.aps.org/doi/10.1103/PhysRevD.84.044030

[15] Grandclément P, Somé C and Gourgoulhon E 2014 *Phys. Rev. D* **90** 024068 URL https://link.aps.org/doi/10.1103/PhysRevD.90.024068

[16] Gralla S E, Holz D E and Wald R M 2019 *Phys. Rev. D* **100** 024018 URL https://link.aps.org/doi/10.1103/PhysRevD.100.024018

[17] Das A, Saha A and Gangopadhyay S 2020 *Eur. Phys. J. C* **80** 180 ISSN 1434-6052 URL https://doi.org/10.1140/epjc/s10052-020-7726-z

[18] Lanza A 1992 *Class. Quantum Gravity* **9** 677-696 URL https://doi.org/10.1088/2F0264-9381%2F9%2F3%2F010

[19] Kleihaus B, Kunz J and Navarro-Lérida F 2002 *Phys. Rev. D* **66** 104001 (Preprint gr-qc/0207042) URL https://link.aps.org/doi/10.1103/PhysRevD.66.104001

[20] Kleihaus B, Kunz J, Mojica S and Radu E 2016 *Phys. Rev. D* **93** 044047 (Preprint 1511.05513) URL https://link.aps.org/doi/10.1103/PhysRevD.93.044047

[21] Delgado J F M, Herdeiro C A R and Radu E 2020 *J. High Energy Phys.* **2020** 180 ISSN 1029-8479 URL https://doi.org/10.1007/JHEP04(2020)180

[22] Villain L, Pons J A, Cerdá-Durán P and Gourgoulhon E 2004 *A&A* **418** 283-294 URL https://doi.org/10.1051/0004-6361:20035619

[23] Abdujabarov A, Atamurotov F, Kucukakca Y, Ahmedov B and Camci U 2013 *Astrophys. Space Sci.* **344** 429-435 ISSN 1572-946X URL https://doi.org/10.1007/s10509-012-1337-6

[24] Cebeci H, Özdemir N and Şentorun S 2016 *Phys. Rev. D* **93** 104031 URL https://link.aps.org/doi/10.1103/PhysRevD.93.104031

[25] Chakraborty C and Bhattacharyya S 2019 *J. Cosmol. Astropart. Phys.* **2019** 034-034 URL https://doi.org/10.1088/2F1475-7516%2F2019%2F05%2F034

[26] Carter B 1969 *J. Math. Phys.* **10** 70-81

[27] Carter B 1973 Black holes equilibrium states *Proceedings, Ecole d'Eté de Physique Théorique: Les Astres Occlus: Les Houches, France, August, 1972* (Les Houches Summer School vol 23) ed DeWitt C and DeWitt B S (Gordon and Breach) pp 57-214

[28] Heusler M 1996 *Black Hole Uniqueness Theorems* Cambridge Lecture Notes in Physics (Cambridge Univ. Pr.)

[29] Gourgoulhon E and Bonazzola S 1993 *Phys. Rev. D* **48** 2635-2652 URL https://link.aps.org/doi/10.1103/PhysRevD.48.2635

[30] Abramowicz M A and Fragile P C 2013 *Living Rev. Relativ.* **16** 1 ISSN 1433-8351 URL https://doi.org/10.12942/lrr-2013-1

[31] Brito R, Cardoso V and Pani P 2015 *Lect. Notes Phys.* ISSN 1616-6361 URL http://dx.doi.org/10.1007/978-3-319-19000-6

[32] Vincent F H, Paumard T, Gourgoulhon E and Perrin G 2011 *Class. Quantum Gravity* **28** 225011 URL https://doi.org/10.1088/2F0264-9381%2F28%2F22%2F225011

[33] Misner C W, Thorne K S and Wheeler J A 1973 *Gravitation* (W. H. Freeman) ISBN 9780716703440. 9780691177793 URL https://press.princeton.edu/books/hardcover/9780691177793/gravitation