Geodesics for the Painleve-Gullstrand form of Lense-Thirring spacetime

Recently, the current authors have formulated and extensively explored a rather novel Painleve-Gullstrand variant of the slow-rotation Lense-Thirring spacetime, a variant which has particularly elegant features -- including unit lapse, intrinsically flat spatial 3-slices, and a separable Klein-Gordon equation (wave operator). This spacetime also possesses a non-trivial Killing tensor, implying separability of the Hamilton-Jacobi equation, the existence of a Carter constant, and complete formal integrability of the the geodesic equations. Herein we investigate the geodesics in some detail, in the general situation demonstrating the occurrence of"ultra-elliptic"integrals. Only in certain special cases can the complete geodesic integrability be explicitly cast in terms of elementary functions. The model is potentially of astrophysical interest both in the asymptotic large-distance limit and as an example of a"black hole mimic", a controlled deformation of the Kerr spacetime that can be contrasted with ongoing astronomical observations.

We emphasize that there is no Birkhoff-like theorem for axisymmetric spacetimes in (3 + 1) dimensions [10][11][12][13][14][15]. The Kerr solution need not, (and typically will not), perfectly model rotating horizonless astrophysical sources such as stars and planets, due to the nontrivial mass multipole moments that such objects typically possess. Instead, the Kerr solution will only model the gravitational field in the asymptotic large-distance regime, a region where the Lense-Thirring spacetime serves as a perfectly valid approximation to Kerr . Given that this variant of the Lense-Thirring metric is a valid approximation for the gravitational fields of rotating stars and planets in the same regime that the Kerr solution is appropriate, there is a compelling physics argument to use the Painlevé-Gullstrand form of Lense-Thirring to model various astrophysically interesting cases [37][38][39][40][41][42][43].
Observationally, apart from its interest in the large-distance asymptotic regime, this Lense-Thirring variant may also be viewed as a "black hole mimic" that can be contrasted with ongoing astronomical observations of various black hole candidates [37,[66][67][68][69][70].
-2 - We note that a competing slow-rotation model has recently been discussed in reference [71]. The trade-off made therein was to improve the integrability properties (the "hidden symmetries") at the cost of sacrificing the global Painlevé-Gullstrand form of the metric.

Killing Tensor and Carter constant
Based on the algorithm presented in two recent papers by Papadopoulos and Kokkotas [72,73], which are in turn based on considerably older results by Benenti and Francaviglia [74], we found the non-trivial Killing tensor [2]: (2.1) Explicitly, we write the metric as [1,2]: (2.2) Then det(g ab ) = −r 4 sin 2 θ, and for the inverse metric we have [1,2]: 3) The (contravariant) non-trivial Killing tensor is [2]: The corresponding covariant form of the Killing tensor, K ab = g ac K cd g db , is then [2]: One can easily explicitly check (e.g., Maple) that ∇ (c K ab) = K (ab;c) = 0.
For any affine parameter λ, the (generalized) Carter constant is now [2]: Without any loss of generality we may choose λ be future-directed, (so dλ/dt > 0). Note that by construction, since it is a sum of squares, C ≥ 0. (For additional recent discussion on general Killing tensors see [75,76].)

Four conserved quantities
In addition to the Carter constant (2.6), in this spacetime geometry we have three other conserved quantities. Two of these come from the time-translation and axial Killing vectors, ξ a = (1; 0, 0, 0) a and ψ a = (0, 0, 0, 1) a , respectively: The final conserved quantity, ǫ, with ǫ ∈ {0, −1} for null and timelike geodesics respectively, comes from the trivial Killing tensor (the metric g ab ): We can greatly simplify these four conserved quantities by rewriting them as: Notice that by construction C ≥ L 2 ≥ 0, and that (dt/dλ) 2 + ǫ ≥ 0.
If ǫ = 0 then, without loss of generality, we can rescale the affine parameter λ to set one of the constants {C, E, L} → 1. It is perhaps most intuitive to set E → 1.
In contrast if ǫ = −1 then λ = τ is the proper time and there is no further freedom to rescale the affine parameter. E then has real physical meaning and the qualitative behaviour is governed by the sign of E 2 + ǫ. Concretely, at least in the case of Carter constant zero, one asks: • Is E < 1? (Bound orbits.) • Is E = 1? (Marginal orbits.) • Or is E > 1? (Unbound orbits).
Forbidden declination range: The form of the Carter constant, equation (3.5), since it is a sum of squares, gives a range of forbidden declination angles for any given, non-zero values of C and L. We require that dθ/dλ be real, and from equation (3.5) this implies the following requirement: Then provided C ≥ L 2 , which is automatic in view of (3.5), we can define a critical angle θ * ∈ [0, π/2] by setting Then the allowed range for θ is the equatorial band: θ ∈ θ * , π − θ * . (3.10) • For L 2 = C we have θ = π/2; the motion is restricted to the equatorial plane.
• For L = 0 with C = 0 the declination is fixed, θ(λ) = θ 0 , and the motion is restricted to a constant declination conical surface. In solving for these unknown functions, one encounters three independent sign choices.

General geodesics
(Note dt/dλ is always taken to be positive.) Consequently, it is useful to define the following quantities (the subsequent physical interpretations are chosen from context for each equation):

Trajectories
For the geodesic trajectories we find the four equations: Here X(r) is the sextic Laurent polynomial In terms of the roots of this polynomial in the generic case we can write In the special case E 2 +ǫ = 0, corresponding to a marginally bound timelike geodesic, the sextic degenerates to a quintic Qualitatively the radial motion can be bounded (if one is trapped between two real roots), or diverge to spatial infinity (if one is trapped above the outermost real root), or be a plunge to r = 0 (if one is trapped below the innermost positive real root). In the immediate vicinity of any real root r i the behaviour will depend on the multiplicity m i of that root. Approximately one has If the root is multiplicity 1, (the generic situation) one has a "bounce" at the turning point (perinegricon or aponegricon -periapsis or apoapsis) at some finite value of the affine parameter: If the root is multiplicity 2, (a somewhat rarer situation), one has an exponential approach to the root: If the root is multiplicity 3 or higher, (an unusual situation; for instance take ǫ → 0 and C → 0), one has a very slow polynomial approach to the root as |λ| → ∞:

Integrating the affine parameter
From equation (4.4), we find This can also be re-expressed as follows: .
Integrals of this type are known as ultra-elliptic integrals [77,78], and date back (at least) to work by Weierstrass and Kovalevskaya in the second half of the 19 th century.
If X(r) were to be cubic or quartic, this would be an ordinary elliptic integral [79]. If X(r) is quintic or sextic, as above, this is an ultra-elliptic integral. More generally, for polynomials of arbitrary order, these would be called hyper-elliptic integrals. Even more generally, these integrals are a sub-class of the so-called Abelian integrals [80]. Generically, equation (4.16) cannot be explicitly integrated in closed form using only elementary functions, hence we cannot analytically invert this relation to find r(λ).
However there is no obstacle in principle to numerical integration to explicitly find the affine parameter λ(r). (Even for the exact Kerr solution one rapidly finds that use of some level of numerical integration is almost unavoidable [81][82][83].) Note that if one is trapped (above or below) any real root of the polynomial X(r) of multiplicity 1, then every time one "bounces" off the root the quantity S r will flip sign, and λ(r) will be double-valued though the inverse function r(λ) will always be single-valued. This is as it should be to guarantee that the affine parameter λ is always continuously increasing with time.
Note that if one is trapped between two real roots of the polynomial X(r), say r min and r max , both of multiplicity 1, then each bounce from r min to r max , or from r max to r min , will advance the affine parameter by some finite amount (the appropriate "period" of the ultra-elliptic integral, now typically called a complete ultra-elliptic integral): Then λ(r) will be multi-valued though the inverse function r(λ) will always be singlevalued. Furthermore, r(λ), while analytically intractable, will at least be known to be periodic, with known periodicity 2 ∆λ. If instead one is approaching a real root of multiplicity 2 or higher, then λ(r) will diverge -one will not actually reach the root for any finite amount of affine parameter lapse -that is r(λ) will asymptote to that higher multiplicity root.
It is straightforward to integrate the first term in the integrand, yielding As before, the remaining integral is an ultra-elliptic integral [77,78], now a different ultra-elliptic integral, and this equation cannot be explicitly integrated in closed form.
Note that if one is trapped between two real roots of the polynomial X(r), say r 1 and r 2 , both of multiplicity 1, and both outside the horizon at r = 2m, then each bounce from r 1 to r 2 will advance the Killing time by some finite amount where now S r = sign(r 2 − r 1 ).
Specifically, if r max = max{r 1 , r 2 } and r min = min{r 1 , r 2 }, then on the upswing from r min to r max one has In contrast on the downswing from r max to r min one has That is The total period is therefore This is again a complete ultra-elliptic integral, now a different complete ultra-elliptic integral.
In this situation t(r) will be multi-valued while the inverse function r(t) will be single -9 -valued. While analytically intractable, r(t) will at least be known to be periodic, with known periodicity 2 T . If instead one is approaching a real root of multiplicity 2 or higher, then t(r) will diverge -one will not actually reach the root for any finite amount of Killing time.

Integrating the declination
As for our equation involving the declination angle θ, first recall the definition of the critical angle θ * as θ * = sin −1 |L|/ √ C . Then from equation (4.6), we find From this we see that is cos −1 cos θ cos θ * = cos −1 cos θ 0 cos θ * + S θ S r √ C r r 0 dr r 2 X(r) . (4.29) Without loss of generality we may allow the geodesic to reach the critical angle θ * at some radius r * , and then use that as our new initial data.
This effectively sets θ 0 = θ * , and gives us the following simplified result: with the last step coming from the fact that cos(· · ·) is an even function of its -10 -argument. That is cos θ = cos θ * cos √ C r r * |dr| r 2 X(r) , (4.31) where the phase is monotone increasing.
One is again reduced to investigating yet another ultra-elliptic integral [77,78], with the declination angle θ oscillating periodically as a function of this phase with period 2π, and with each "bounce" from r min to r max advancing the phase of the cosine by an amount ∆(phase) = √ C rmax r min dr r 2 X(r) . (4.33) As before, this equation cannot be explicitly integrated in closed form.

Integrating the azimuth
We now finally consider the ODE for the evolution of the azimuthal angle: dφ/dλ. (This particular sub-case is considerably messier than the previous ones.) Using equations (4.7) and (4.4) we find (4.34) Consequently where as per our previous discussion we have S r = sign(r − r 0 ). The last term appearing here is again an ultra-elliptic integral [77,78]. In contrast, the penultimate term is somewhat worse, because the integrand contains θ(r), the overall integral is an iteration of an ultra-elliptic integral.
-11 - We can explicitly integrate the first integral in closed form: Thence, in general As before, this equation being ultra-elliptic means it cannot be explicitly integrated in closed fully analytic form, at least not in terms of elementary functions.

Summary of generic geodesic evolution
Overall, while these generic geodesic equations cannot be integrated in closed and fully analytic form, they are still integrable in the formal technical sense. In terms of complete ultra-elliptic integrals some of the key quantities of interest (both mathematically and physically) are the "periods": Typically, the "periods" of these ultra-elliptic integrals will be incommensurate.
If further specific constraints are now imposed, then these equations can indeed be integrated in closed form. In the next section we explicate both the null and timelike geodesics for when C = 0 in closed fully analytic form. We are able to recover the "rain" geodesics from [1], as well as present a simple derivation of both the "drip" and "hail" geodesics.

Geodesics with Carter constant zero
If the Carter constant is zero then the geodesic equations simplify radically. Firstly, if C = 0 then from the manifest positivity of we see that we must have both L = 0 and dθ/dλ ≡ 0. The condition that L = 0 physically constrains the geodesics to the trajectories of ZAMOs (zero angular momentum observers), whilst dθ/dλ = 0 =⇒ θ(r) = θ 0 ; some constant θ 0 .
Furthermore our expression for X(r) significantly simplifies since now: We find that the four trajectory equations reduce to: The form of these equations suggests it would be particularly useful to separate our analysis of null geodesics (photons) from timelike geodesics (massive particles).

Null geodesics (photons) with Carter constant zero
For null geodesics (photons) with Carter constant zero, we have the following conditions: Furthermore, without loss of generality, we can rescale the affine parameter to set E → 1, so that X(r) → 1. From equation (4.16), we now find That is, we find the simple linear relation Thus in this particular situation r is an affine parameter. Ingoing geodesics will crash into the central singularity in finite affine time, whereas outgoing geodesics can emerge from the horizon (r = 2m) at finite affine time (which we shall soon see corresponds to minus infinity in Killing time). The apparent asymmetry between ingoing and outgoing null geodesics is a side-effect of the initial choices made in setting up the Painlevé-Gullstrand coordinate system. (Did one choose an ingoing or outgoing Painlevé-Gullstrand coordinate system?) Furthermore, equation (4.20) for t(r) now reduces to This integrates explicitly to This can also be rewritten as Which form one uses is really a matter of taste.
Finally, from equation (4.37) we also have This integrates explicitly to Notice that while these equations are rather complicated, they are all fully explicit and given in terms of elementary functions. We also note that these equations have sensible limiting behaviour; as r → r 0 , we have that λ(r), t(r), φ(r) → λ 0 , t 0 , φ 0 respectively.

Timelike geodesics (massive particles) with Carter constant zero
For timelike geodesics (massive particles) with Carter constant zero, we have the following conditions: The form of the polynomial X(r) suggests that we should split our analysis into three cases. For the first case we set E = 1 (since in this case X(r) reduces significantly; X(r) → 2m/r), for the second case we set E > 1, and for the third case we set E < 1. Physically, geodesics with E = 1 represent particles that have zero radial velocity at spatial infinity (these are marginally bound geodesics). Geodesics with E > 1 represent particles that have non-zero velocity at spatial infinity (these are unbound -15 -geodesics). Geodesics with E < 1 represent particles that are in bound orbits, and so never escape to spatial infinity. For ingoing geodesics, these correspond to the "rain" (E = 1), "hail" (E > 1), and "drip" (E < 1) geodesics.

Marginal geodesics E = 1
Setting E = 1, our conditions reduce to: So, for our expression for λ(r) we find Our general expression (4.20) for t(r) reduces to We may now compute the somewhat unwieldy integral r r 0 giving the following explicit form for t(r) It is worth noting that for S r = −1 (corresponding to ingoing geodesics) we find the particularly simple result whilst for S r = +1 (corresponding to outgoing geodesics) we obtain The apparent asymmetry between ingoing and outgoing timelike geodesics is a sideeffect of the choices made in setting up the Painlevé-Gullstrand coordinate system. Note that in this situation the formulae for t(r) are J-independent, so this apparent asymmetry is already present for Schwarzschild spacetime in Painlevé-Gullstrand coordinates. The usual choices allow ingoing geodesics to penetrate the horizon and crash into the central singularity in finite Killing time, while outgoing geodesics need an infinite amount of Killing time to escape from the horizon at r = 2m and so get 'stuck'.
-17 -Our general expression (4.37) for φ(r) reduces to We may integrate the remaining integral in closed form This gives the fully explicit form for φ(r) as: ) .

(5.26)
It is now straightforward to recover the "rain" geodesics explored in [1], which model a ZAMO dropped from spatial infinity with zero initial radial velocity. This physical scenario requires C = 0, L = 0, ǫ = −1, and E = 1, as above, and furthermore is an ingoing retrograde geodesic (corresponding mathematically to fixing S r = S φ = −1). These geodesics are retrograde since we must give them an initial nonzero angular velocity in the retrograde direction at spatial infinity in order for L = 0 to hold along the length of the geodesic. From equations (5.22) and (5.26), we find the explicit "rain" geodesics: If we differentiate with respect to r, we find which are exactly the equations defining the rain geodesics as given in reference [1].

Unbound geodesics E > 1
For E > 1, let us first restore the full list of conditions from equation (5.15): Notice that for E > 1, r is a priori unconstrained; r ∈ (0, +∞).
From our general result (4.16) for λ(r) we have Conducting a Taylor series expansion around E = 1, we find In the limit where E → 1, we see that our expression for λ(r) simplifies to (5.17), as expected.
-19 - The unwieldy integral in our general expression (4.20) for t(r) now becomes r r 0 This yields the following fully explicit result for t(r) Conducting a Taylor series expansion around E = 1 gives (5.37) Hence in the limit E → 1, we have which is just equation (5.21), as expected.
Lastly, the somewhat unwieldy integral in our general expression (4.35) for φ(r) is now given by so we find that φ(r) reduces to If we now make the direct substitution E = 1, we find which is just equation (5.26), as expected.
Overall, while these equations of motion are rather long, they are fully explicit and have appropriate limits when E → 1.
In all generality, geodesics with E > 1 are unbound geodesics. When (as herein) the Carter constant is zero, the ingoing geodesics are the so-called "hail" geodesics, modelling a ZAMO fired in from spatial infinity with a nonzero initial velocity. These geodesics are ingoing retrograde, and are hence given by and

Bound geodesics E < 1
Once again we repeat the full list of conditions from equation (5.15): If we let E < 1, then X(r) = 0 has a unique root at r = 2m 1−E 2 , and in order to keep X(r) real, we see that we have the constraint that r ∈ (0, 2m 1−E 2 ]. In particular r ≤ 2m/(1 − E 2 ). Physically, geodesics with E < 1 are gravitationally bound (and, because C = L = 0, must eventually crash into the central singularity). These particular bound geodesics are the so-called "drip" geodesics, corresponding to a ZAMO being dropped from some finite r * with zero initial velocity.
At intermediate stages of the computation it is useful to use the identity ln(x + iy) = 1 2 ln x 2 + y 2 + i cos −1 x For S r = +1 these drip geodesics will crash into the central singularity r = 0 in finite affine time (5.48) These "drip" geodesics are qualitatively (not quantitatively) somewhat similar to the "rain" geodesics given by equations (5.27) and (5.29).

(5.51)
Ingoing geodesics S r = −1 (in this context the "drip" geodesics) will crash into the central singularity in finite Killing time There are of course many other ways of rearranging this result.
Finally for φ(r) the result (5.40) can be extended to the region E < 1 without alteration.

Conclusions
From this discussion we have seen that, once given the non-trivial Killing tensor for the Lense-Thirring spacetime, we can extract the Carter constant; the fourth constant of the motion. Then the geodesic equations become integrable at least in principle (in terms of ultra-elliptic integrals). This allows us to formally solve for myriads of general geodesics. However, we saw that in full generality, we could not explicitly integrate the equations of motion in closed form, the ultra-elliptic integrals are not "elementary". Only when imposing further conditions, such as Carter constant zero, could we then explicitly integrate the equations of motion in algebraically closed form.
-24 - The explicit geodesics given in this discussion are quite tractable and can be applied to a number of astrophysically interesting cases. For example, the Carter constant zero geodesics with E < 1 are the "drip" geodesics of the spacetime, while the Carter constant zero geodesics with E = 1 are the "rain" geodesics, and the Carter constant zero geodesics with E > 1 are the "hail" geodesics of the spacetime.