Kerr-fully Diving into the Abyss: Analytic Solutions to Plunging Geodesics in Kerr

We present closed-form solutions for plunging geodesics in the extended Kerr spacetime using Boyer-Lindquist coordinates. Our solutions directly solve for the dynamics of generic timelike plunges, we also specialise to the case of test particles plunging from a precessing innermost stable circular orbit (ISSO). We find these solutions in the form of elementary and Jacobi elliptic functions parameterized by Mino time. In particular, we demonstrate that solutions for the ISSO case can be determined almost entirely in terms of elementary functions, depending only on the spin parameter of the black hole and the radius of the ISSO. This extends recent work on the case of equatorial plunges from the innermost stable circular orbit. Furthermore, we introduce a new equation that characterizes the radial inflow from the ISSO to the horizon, taking into account the inclination. For ease of application, our results have been implemented in the KerrGeodesics package in the Black Hole Perturbation Toolkit.


Introduction
Analysing the geodesic structure of a spacetime is one of the foundational means of understanding it in its unperturbed state.This is apparent in the fact that the geodesics of the Kerr spacetime have been extensively studied since its original derivation [1].A critical step in the calculation of the geodesics in Kerr was the discovery of a fourth constant of motion, the Carter constant Q [2], which along with the mass shell condition, the conserved energy E and the angular momentum L allow for the full integrability of the geodesic equations of motion.While exact solutions for various special cases had been derived (e.g, see [3] or [4]) no real effort was made to tackle the generic case until the start of the 21st century [5].Without getting explicit solutions for the geodesics themselves, [6] found exact expressions for the orbital frequencies with respect to coordinate time of generic bound orbits.The introduction of the Mino time parameterisation [7] subsequently allowed for the full decoupling of the problem in a much simpler manner to previous approaches [2].This opened the door for finding analytic solutions to generic bound geodesics in Kerr [8] as a system of piecewise smooth functions, which were then notably simplified through analytic continuation [9].
The full class of analytic solutions to Kerr-de Sitter and Kerr-anti-de Sitter spacetimes has been found generically [10].However these solutions are presented in the form of Weierstrass elliptic and hyperelliptic Kleinian functions, which can make them cumbersome to deal with.This work also provided part of the motivation for deriving a more explicit solution for null geodesics in the exterior of Kerr in terms of Jacobi elliptic functions [11].Recent work has subsequently derived these equations for the specified case of bound null geodesics [12], relevant in pursuit of forming tight constraints on Black Hole parameters by observations of the photon ring using the Event Horizon Telescope [13].
Recently, [14] has derived a much simplified analytic solution for the special case of equatorial plunging timelike geodesics which asymptote from the innermost stable circular orbit (ISCO).In this work, they also provide a simple expression for the equatorial radial inflow from the ISCO relevant to the study of accretion disk dynamics.In practice, these systems will not always be confined to the equatorial plane, motivating the generalisation of this result to the inclined (i.e.non-spin aligned or precessing) case, as we will do in this paper.In the interest of completeness we also provide solutions for plunges starting from a general (inclined and eccentric) last stable orbit (LSO), i.e. asymptoting from a generic unstable spherical orbit (USO).
Our work on generic plunges is further motivated by the fundamental role that geodesics play in the gravitational self-force approach to solving the relativistic dynamics of binary black holes [15].In this approach the dynamics are expanded in powers of the mass-ratio between to two black holes.In this scheme, the zeroth approximation to the motion of the lighter secondary component is given by a geodesic in the Kerr geometry generated by the (heavier) primary black hole.At higher orders, this motion is corrected by an effective force term, the gravitational self-force, causing the system to evolve.During the inspiral phase this evolution can be solved using a 2-timescale formalism [16][17][18], adiabatically evolving the system along a sequence of bound geodesics.The 2-timescale formalism breaks down as the system approaches the LSO around the black hole, where it enters a transition regime governed by a new intermediate timescale [19][20][21][22][23].In the asymptotic regime beyond the transition, the motion is again well described by a perturbed geodesic, now of the plunging variety.The timescale for the plunge is much shorter than radiation reaction timescale associated with the self-force.Consequently, the dynamics in this regime are dominated by the geodesic term.In practice, full inspiral-transition-plunge trajectories are formed by asymptotically matching an inspiral trajectory and a plunging geodesic to the jump obtained from analysing the transition regime [20,24].
Development of the gravitational self-force formalism was originally motivated by the need for producing accurate gravitational waveforms for the observation of extreme mass-ratio inspirals (EMRIs) with the planned space-based gravitational wave observatory LISA.EMRIs are expected to have mass-ratios of order 10 −5 , and will therefore spend hundreds of thousands of orbits in the strong field regime of the inspiral phase.Consequently, the handful of orbits represented by the transition and plunge phases are generally expected to provide a negligible contribution to the total signal.However, over time, evidence has mounted [25][26][27][28][29][30][31][32][33][34] suggesting that the self-force formalism can produce accurate results at much higher mass-ratios, and possibly even for comparable mass binaries.In these regimes the plunge and associated ringdown form a much more significant portion of the waveform.Consequently, the realisation that waveforms from the self-force formalism may be usable in the more comparable mass regime has led to a renewed interest in modelling the transition and plunge phases [24,[35][36][37][38][39][40] and the gravitational waves produced during the plunge and subsequent ringdown [41][42][43].So far this effort has focused mostly on special cases involving quasi-circular (possibly precessing) inspirals.As the work progresses towards fully generic inspirals, there is a need for generic solutions for the plunge geodesics.This paper provides the latter by solving for generic plunging geodesics in an easy to evaluate form.
The layout of this paper is as follows, in section 2 we begin with an introduction on the geodesic equations of Kerr and provide explicit definitions for the related conserved quantities.In section 3 we then focus our attention on the special case of plunging timelike geodesics which asymptote from the innermost stable precessing circular orbit (ISSO).‡ From these equations we determine the exact and two approximate expressions for the rate of radial inflow from the ISSO to the horizon.We then go on to determine the fully analytic solutions to these geodesic equations for ISSO plunges.In section 4 of this work, we determine novel, fully analytic, expressions for generic timelike plunging geodesics in the Kerr spacetime, in terms of elementary and (Jacobi) elliptic functions.These generic plunges are presented in a manifestly real form such to be easily implementable, supporting current work in the self-force community on the aforementioned transition to plunge.In relation to solutions of special classes of geodesics we also provide the solutions for plunges asymptoting from an USO in Appendix A. Finally, we have implemented these solutions in the KerrGeodesics package of the Black Hole Perturbation Toolkit.We work in geometric units where G = c = 1.

Geodesic Equations
We work in modified Boyer-Lindquist coordinates (z = cos(θ)) and denote the mass and spin of the black hole as M = 1 and a, respectively.We further use the standard definitions, Σ = (r 2 + a 2 z 2 ) and ∆ = (r − r + )(r − r − ).The inner and outer event horizons are given by (1) ‡ In this paper, we use the acronym ISSO to refer to the general case of the last stable circular orbit for spherical (i.e.inclined precessing) orbits.The term ISCO will be used exclusively for the special case of circular equatorial last stable orbits.
The symmetries of the Kerr geometry provide two constants of motion, the conserved energy and angular momentum, which are given as respectively.Here g µν is the Kerr metric tensor in modified Boyer-Lindquist coordinates.
There also exists a third constant of motion Q known as the Carter constant [2] which arises in some families of spacetimes exhibiting Type D symmetry [44].These conserved quantities arise as a result of the existence of certain Killing tensors, which are defined to satisfy the condition ∇ (a K bc) = 0.In Kerr, a tensor satisfying this conditions can be found, and is given by, Here l µ and n ν are the principal null vectors of the Kinnersly tetrad [45] defined by and . The Carter constant is then given by, ( The Kerr metric also satisfies the Killing tensor condition, giving rise to a final constant of motion, g µν u µ u ν = −1, also known as the mass shell condition,.Taking these four conserved quantities the equations of geodesic motion in Kerr are given by [2], dz dλ where we have specialised to working in Mino(-Carter) [7] time defined by By taking advantage of the Mino time parameterisation the equations of motion completely decouple and can be solved hierarchically.This is done by first solving for r(λ) and z(λ) then naturally solving the equations in the form t(r, z, λ) = t r (r) + t z (z) − aEλ and ϕ(r, z, λ) = ϕ r (r) + ϕ z (z) + aLλ.The solutions for r(λ) and z(λ) can be directly substituted to obtain t(λ) and ϕ(λ).
Analysing the radial equation in its explicit form (6), one can see that the distinction between bound and plunging orbits is fully determined by the root structure of the fourth order polynomial R(r).In particular, plunges occur when we have bound motion between some roots r i , r j of R(r) with r i < r − < r + < r j .In this work, we restrict to geodesics with E < 1, which implies that the radial potential R(r) is negative in the limit r → ∞, ensuring that the geodesic is bound to the Kerr black hole.Moreover, since R(r ± ) > 0 it guarantees that there is at least one real root of R(r) outside r + .A less obvious implication comes from the polar equation (7).Rewriting the polar potential as it becomes apparent that when E < 1 the polar equation has solutions with −1 < z < 1 if and only if Q ≥ 0. For the radial potential this implies that R(0) ≤ 0, and consequently that there exists at least one real root of the radial equation between r = 0 and the inner horizon.We thus find that for any values of E < 1 and Q ≥ 0, there exists a plunging geodesic.It is important to note that the case of E > 1 with Q < 0 also contributes to a small subset of parameter space for which a real solution is allowed.Its dynamics are quite interesting as by a brief analysis of the equations of motion one can see that this would give a test particle which plunges from infinity, poloidally oscillating around some z 0 ̸ = 0. Generically, a plunging geodesic will eternally oscillate between two turning points of the radial potential and return to the same radial point after a finite amount of Mino (and proper) time.Geometrically, this corresponds to a geodesic diving into the Kerr black hole and passing through the two horizons before being scattered back out, passing the horizons in reversed order and exiting in a different asymptotically flat region of the maximally extended Kerr solution, as shown in the Penrose diagrams in Fig. 1.
Exceptions to this picture occur for equatorial trajectories and when one of the bounding roots has a multiplicity greater than one.In the first case Q = 0, one finds that r = 0 is a root of the radial equation.If this is the inner turning point of the geodesic, the plunge will end on the singularity after a finite amount of Mino (and proper) time.In the second case, approaching a root with higher multiplicity takes an infinite amount of Mino (and proper) time.In Section 3, we will see the physically relevant case where the outer turning point of the solution is a triple root, i.e. lies on the ISSO, where the radial solution takes a particularly simple form generalising the result of [14].The right panel in Fig. 1 shows this edge case trajectory in a Penrose diagram.In Appendix A, we also treat the special case where the outer turning point is a double root, generalising some of the results in [46].
In the generic case, we know at least two of the four roots of R(r) are real.The other two roots are either both real or both complex.If the other two roots are real, they come in a pair that lies either entirely outside the outer turning point, entirely between the inner turning point and r = 0, or entirely in the r < 0 region.In the first of these cases, the plunging orbit exists inside of a normal bound orbit, a case sometimes referred to as a "deeply bound" orbit.The solutions for cases with 4 real roots turn out to be a straightforward generalisation from the solutions of [8,9].The derivation of the solution for the complex case will turn out to be substantially more involved.

The Innermost Precessing Stable Circular Orbit
Before determining the solutions to fully generic plunges in Kerr we begin by solving for plunges which asymptote to the ISSO.In this case the radial potential R(r) is imbued with a triple root at the ISSO radius r = r I .Two of these roots come from the fact that the ISSO must be a (precessing) circular orbit and the additional root arises from the fact we are looking specifically at the innermost of these orbits meaning R(r) must also inflect at this point.As a result we determine a radial equation of the form By equating Eq. ( 6) with Eq. ( 12) one immediately obtains the result, The roots of the polar equation can be found to be given by [11], Finally, we define as a quantity which will recurrently show up throughout this work.§

Determining the Conserved Quantities
Naturally geodesics asymptoting to the ISSO must also share the same constants of motion as the ISSO.We can therefore identify a plunging geodesic of this type with those of a particular ISSO, which has two degrees of freedom.These two degrees of freedom can be set by picking a black hole spin (a) and maximum orbital angle of inclination ).This then determines a unique r I .One can also invert this relation to set the parameterisation in terms of (a, r I ).Making this choice one finds for each value of a there exists a range of allowed r I 's each of which corresponds to a unique inclination either in prograde or retrograde.The innermost and outermost of these quantities correspond to the equatorial ISCOs in prograde and retrograde respectively.
Forms of these bounds have been known in the literature for some time [47].If one wishes to parameterise by inclination, one can use the KerrGeodesics package in the Black hole perturbation theory toolkit [48] to find E, L and Q parameterised by (a, θ inc ), where θ inc runs from 0 for equatorial prograde orbits to π for equatorial retrograde orbits [49].
In parameterising the conserved quantities by (a, r I ) we begin with the expression for the marginally stable spherical orbits Q written in terms of r I [50] which gives, § Note that the definition of k z (and later k r ) differs from the conventions used in [9].
Equating the remaining coefficients between Eq. ( 6) and Eq. ( 12) we obtain the equations, Where the ± is determined by whether or not the r I picked corresponds to a prograde or retrograde orbit respectively.The correct sign is determined by the condition, where, defining κ = √ a 2 − 2r + r 2 , the value of L root is given as the root of the function which is real and closest to r = 6.Eq. ( 21) has been given in a form such that the root can be found numerically to high precision which is not the case for Eq.(19).It is worth noting that as E, L, and Q have all been determined in terms of r I , and z 1 is the root that defines the maximum range of oscillation allowed for a given r I , our solution immediately defines a spacelike surface (r, z, ϕ) = (r I , z 1 (r I ), u) for r I ∈ (r Imin , L root ) and u ∈ (0, 2π) inside which no stable spherical orbits can exist.Taking our exact solution for this spacelike surface to the extremal limit also provides us with the ISSO surfaces found in the near-horizon extremal Kerr geometry previously found in both [51] and [52].

Radial Inflow
We are now ready to begin analysing the physical consequences of extending the results of [14], regarding equatorial ISCO flow to inclined orbits.The exact solution to this equation includes functional dependence on certain Jacobi elliptic functions.In order to simplify the results we provide two approximate forms of the inflow.The first, approximated to a modified form of the equatorial inflow equation which removes all dependence on Jacobi elliptic functions.The second, a polar averaged form which simplifies the functional dependence on the Jacobi elliptic functions to a constant dependence for any given set of parameter values.We find the modified equatorial flow to be given by, Plot of differing radial inflow equations from the ISSO to horizon, for parameter values (a, r I ) =(0.9,4) for plot range (r + , r I ) on the radial axis.The lightly oscillating purple line corresponding to the exact solution exhibiting small oscillations over each polar period.The blue and orange dashed lines correspond to the equatorial and polar averaged approximations respectively.The Radial flow measured on the y axis is given by the dimensionless quantity dr/dt for each of the three flow equations.
In the equatorial limit (Q = 0), this reproduces the result of [14], but gives an improved representation of the radial inflow for inclined disks where Q, E ,and L are given by their true values Eqs. ( 17),( 18) and ( 19) respectively.Next we go on to determine the form of the radial inflow when averaged over the polar period, this is found to be Where K(•) and E(•) are the complete elliptic functions of the first and second kind respectively.We have defined the polar average as follows, given some function f (λ) that (partially) depends on λ through z such that f (λ) = F (r(λ), z(λ), λ), we take where is the polar period.For example when we have an equation of the form t(r, z, λ) = t r (r) + t z (z) + aLλ then ⟨t⟩ z (λ) = t r (r(λ)) + ⟨t z ⟩ z (λ) + aLλ.Here the purpose of taking this average is to integrate out the oscillatory dependence and isolate the secular dependence in λ of terms originally containing polar dependence.Finally, we give the exact equation for the radial inflow as, Where sn(•|•) is the Jacobi sine function.In Fig. 2, which depicts the radial flow from the differing equations, we see an improvement of accuracy in the polar averaged approximation over the equatorial approximation.We also provide a brief error analysis comparing our approximated solutions to the true solution Eq. ( 26) over the entire parameter space.Here we parameterise our space using (x inc , a) where x inc = cos θ inc such that (x inc , a) ∈ ({−1, 1}, {0, 1}).From Fig. 3 we see that over the parameter space the equatorial approximation confines the error to be below 5% whilst the polar averaged form provides a bound of 2% maximum relative error.In reality this is a generous upper bound and for the majority of parameter space the error of the approximated particle inflow will be notably smaller, as can be seen in Fig. 3.

Solutions to the Equations of Motion
At this stage we are now ready to solve the equations of motion.In the generic case we expect the solutions to arise in terms of elliptic functions.In the ISSO case however, we find that the triple root significantly simplifies the terms with radial dependence to the form of elementary functions.For convenience the solutions for the entirety of the ISSO case are given with initial conditions (t(λ), r(λ), z(λ), ϕ(λ))| λ=0 = (0, r 4 , 0, 0).The full solution can be easily reconstructed using Eq. ( 27).

Polar Equation
The polar solution is unchanged from the generic bound case [9] and given by z(λ) = z 1 sin ξ z (λ).
Here we have defined, where am(•|•) is the Jacobi amplitude function.Conveniently for the remaining equations of motion, the polar and radial dependence in the time and azimuthal equation fully decouple from each other in Mino time.This means we are only required to resolve the radial component as the polar part will simply remain the same as has already been found for bound orbits in [8,9].The polar dependence still remains in the form of elliptic integrals where F(•), E(•) and Π(•) are the elliptic integral of the first, second and third kind respectively.

Azimuthal Equation
Going on to solve the azimuthal component we first rewrite Eq. ( 9) as We then integrate each of the terms individually.The component comprising of r dependence is found to be given by, where the arrow notation denotes taking the other term within the shared brackets and swapping all occurrences of r − with r + .The component with z dependence is then given by, The polar average form of this solution can also be found and dramatically simplifies the functional dependence, where Π(•, •) is the complete elliptic integral of the third kind.Thus, we find the full azimuthal solution to be given by

Time Equation
We perform a similar separation in the integration of the equation of motion for coordinate time giving, Next we find the polar dependent part of the time solution to be given by, and the polar averaged form of this solution is given by, This is the expression that was necessary to derive the polar averaged radial flow in Eq. (23).Putting this all together we then find the full time solution to be given by We find that all of the novel radial integrals which are required to be solved can be done so without too much issue through the use of partial fractions.They can then be fully analytically continued by applying some simple trigonometric substitutions.
Additionally one can follow this approach to find the solution for proper time as a function of Mino time along a plunging ISSO geodesic, From Fig. 4 it can now be seen explicitly that enforcing the triple root places us in a regime where these geodesics asymptote from the ISSO and subsequently plunge in through the horizon.This showcases a number of properties of the plunge in Boyer-Lindquist co-ordinates.The structure of these geodesics is more easily depicted when orthogonally projected onto azimuthally co-precessing and poloidally co-rotating planes fixed to a particle as it follows the geodesic.In Fig. 5 we see that the azimuthal coordinate diverges at both the inner and outer horizons.Analysis of the coordinate time solution shows it also diverges at these points, this occurs for a similar reason as to the coordinate time divergence at the horizon in Schwarzschild coordinates for a non- Co-rotating orbital plots of plunging geodesics which asymptote to the ISSO in the infinite past with (a, r I ) = (0.9, 2.6) .Left Image: orthogonal projection of ISSO plunge onto the co-precessing azimuthal plane.Right Image: orthogonal projection of ISSO plunge onto co-rotating polar plane.The cyan and orange frames represent the projection onto the planes seen in 4 as they track a particle following the geodesic trajectory.
spinning black hole, i.e. due to the infinite redshift.In much the same way as we have an infinite redshift on the horizon one can intuitively think of this azimuthal divergence occurring as a consequence of the infinite redshift on the horizon forcing the geodesics to also co-rotate with the black hole an infinite number of times before passing through the horizon.Naturally this is only a co-ordinate singularity and the process occurs in finite Mino and proper time.Due to this discontinuity care must be taken in the choice of branch between r + and r − .The solutions depicted in 4 and 5 seem to invert between the horizons but it can be checked that our solutions conserve E, L and Q throughout their parameterisation, confirming we have chosen the correct branch.In the equatorial (Q = 0) limit, these solutions are found to agree with the results of [14].

Fully Generic Plunging Orbits
We now move on to the case of generic timelike plunging geodesics in Kerr with no restriction on the root structure of the effective radial potential.In this regime the integrals we need to compute are notably more involved, separating into two primary classes.The first being when the radial potential admits four real roots and the second being when the radial potential admits two real and two complex roots.In the case with R(r) having four real roots (r 4 < r 3 < r 2 < r 1 ) with r 4 < r − < r + < r 3 < r 2 < r 1 the solution can be directly borrowed from the bound orbit case [8,9] with the substitution r 1 ⇐⇒ r 3 and r 2 ⇐⇒ r 4 .In addition, in the case of 4 real roots (r 4 < r 3 < r 2 < r 1 ) with r 4 < r 3 < r 2 < r − < r + < r 1 then the solutions can again be found from [8,9], this time with no modification.Naively, in the case of two real and two complex roots (the interesting case for self-force plunges) one may think it to be sufficient to simply take the previously found bound solutions and analytically continue them to allow for complex values of the radial roots.With care this can be done to give correct answers, however intermediate terms in the evaluations give rise to large cancellations between complex numbers effecting numerical accuracy and evaluation speed.To find these solutions in a manifestly real, easy to evaluate, and practical manner, we are required to begin the procedure of solving the complicated elliptic integrals from scratch.

Generic Plunge Orbits with Two Complex Roots
In calculating the integrals for this case we follow the procedure outlined in [53] and give an overview of the steps involved.We begin by again noting that only the radial components of each of the equations need to be solved as the other remaining terms are identical to those already found in the ISSO case above.Recall we now have the general expression R(r) = (E(r 2 + a 2 ) − aL) 2 − ∆(r 2 + (aE − L) 2 + Q) where E, L and Q are all independent quantities.At this stage the integrals of concern are given by Although these integrals look to be of the same form as those we just solved for in the ISSO case without much mention, the cause for the substantial increase in complexity is due to the fact that R(r), in general, no longer contains any double or triple roots.This forces us to much more carefully consider the elliptic integrals at hand.The key idea in the procedure we wish to apply is to reduce the integrals of Eq. 43 such that we must only solve integrals of the form, We solve these integrals in terms of the radial co-ordinate r, which can then be parameterised in terms of Mino time by inverting the solution for λ(r).We begin solving these integrals by first continually applying partial fractions to the radial parts of the ϕ and t equations until arriving at the forms, At this point we can now concentrate on calculating the four elliptic integrals defined in Eq. (44).We do this my applying a transformation given for the case of two complex roots in [53].This procedure begins by letting R(r) have two real roots r 2 < r 1 and two complex roots r 3 and r 4 .We then rewrite R(r) in the form ) where ρ r = ℜ(r 3 ) and ρ i ℑ(r 4 ).Further we define , and Next, motivated by the tables provided in [53], we make the substitution in the integrals Eq. ( 44) of the form, Applying this transformation to Eq. ( 44) and again repeatedly applying partial fractions we find each integral reduces to a sum over elliptic integrands and rational polynomials.The solutions we then obtain are only analytic on the range r ∈ (r 2 , r 1 A+r 2 B A+B ) where for convenience we have set the initial conditions to be given by (t(λ), r(λ), z(λ), ϕ(λ))| λ=0 = (0, r 2 , 0, 0, ).Next we analytically extend these solutions through the point r = r 1 A+r 2 B A+B by use of trigonometric and elliptic substitutions providing a fully analytic solution on the range r ∈ (r 2 , r 1 ).The fully analytical solution to these integrals is then given by dr R(r) = 1 and, Where we have defined, , and In the above we have suppressed the explicit dependence of ξ r on λ in the interest of readability.The arctan function seen in Eq. ( 51) is the two argument arctan function which tracks the sectors of the numerator and denominator.

Solutions to Equations of Motion
The solution for Eqs. ( 49)-( 52) for the elliptic integrals can then be readily substituted back into Eq.( 45) and Eq. ( 46) giving the full form of the solutions to generic plunges.We present all solutions parameterised in terms on Mino time (λ), which is done by inverting the solution to Eq. ( 49) to obtain r(λ) then substituting the solution for r(λ) everywhere r appears in Eqs. ( 50)- (52) .These solutions are provided in a fully analytic form.The radial equation is first found by inverting the solution for Eq. ( 49) to give, The solutions to the polar equation remain the same as for the ISSO case.Taking Eq. ( 45), the solution to the azimuthal equations of motion can then immediately be found from the solutions of Eq. ( 44) Similarly, from Eq. ( 46) we can now obtain where both ϕ z and t z can be taken from the ISSO case.The solution for proper time as a function of Mino time can also be found as, τ r (λ) =I r 2 (λ), and Having obtained the full set of solutions for generic plunges we plot the spatial component depicting the orbital evolution of the generic plunging geodesics Fig. 6.More informatively, the orthogonal projection onto the azimuthally co-precessing and poloidally co-rotating planes in Fig. 7 again show the divergences at either horizon in the ϕ coordinate.Importantly, once we have constructed these solution, we check that E and L are indeed still conserved by evaluating −u ν g µν ∂ ∂t ν and u ν g µν ∂ ∂ϕ ν explicitly.We confirm this is the case for all values of λ, not only acting as a consistency check of our equations, but also showing we have selected the correct branches of the solution between each of the horizon divergences.As a final consistency check, we substitute our solutions back into the equations of motion for both the ISSO and generic case and find that our solutions do in fact solve the original equations.Finally, from our solutions the radial and polar frequencies with respect to Mino time for a generic plunge can be obtained, , and

Discussion
In this work we have obtained the closed form analytic solutions for generic bound plunging geodesics in Kerr.We have also paid particular attention to the special edge case of plunges starting asymptotically from the innermost stable precessing circular orbit, in which the solutions take a particularly simple form.This generalises the result of [14] to inclined motion.The general expression for the inflow in this case, is more involved due to oscillations coming from the polar motion.Therefore we have provided two simplified approximations for the inflow rate.We have also found that the geodesics asymptoting from the ISSO can be parameterised purely in terms of black hole spin and the radius of the ISSO.We expect these solutions to have applications in the modelling of accretion flow in Kerr spacetimes.In addition to the special case of ISSO plunges, we expect the provided generic solutions for plunging geodesics to find practical use in modelling the inspiral of binary black holes, since in the small mass-ratio limit they will describe the final phase of the inspiral before merger.As small massratio methods are applied to more equal mass systems, including this phase becomes increasingly important.For ease of application we have incorporated our results in the KerrGeodesics package in the Black hole perturbation toolkit [48].
In this work we have restricted attention to "bound" plunging geodesics, i.e. geodesics with E < 1.In principle, the explicit solutions given here can easily be extended to the case of "direct plunges" coming from infinity and falling into the black hole, or deeply bound plunges inside a scattering geodesic.In both cases, this involves taking the expressions in this paper in terms of the outermost real root and analytically continuing it past positive infinity to negative values (e.g. by considering the reciprocal root).However, this will not give all E > 1 plunging geodesics, which also includes so-called "vortical" geodesics with negative values of Q [51] leading to qualitatively different behaviour of the polar solutions.

Figure 1 .
Figure 1.Penrose diagrams of the maximally analytic extension of the Kerr metricshowing different plunge trajectories.Left Image: A generic plunging orbit whereby the solutions oscillate between the largest root of R(r) which is less than r − and the smallest root of R(r) greater than r + .These trajectories correspond to geodesic motion entering the black hole, reaching a turning a point, and exiting the white hole into a different asymptotically flat region.Each point on the dashed lines in region I and I ′ correspond to a spacelike 2-surface of constant radius corresponding to the outer turning point.Right Image: A plunge which asymptotes to r I as λ → ±∞, where γ(λ) represents the geodesic in question.On the Penrose diagram these trajectories begin at past spacial infinity i − and end at future spacial infinity i + in a different asymptotically flat region.Each point on the dashed lines in region I and I ′ correspond to a spacelike 2-surface of constant radius corresponding to the ISSO radius r I

1 3
and B = 3 + a 2 + 9−10a 2 +a 4 A + A the range of possible r I values for a given a are,

Figure 2 .
Figure2.Plot of differing radial inflow equations from the ISSO to horizon, for parameter values (a, r I ) =(0.9,4) for plot range (r + , r I ) on the radial axis.The lightly oscillating purple line corresponding to the exact solution exhibiting small oscillations over each polar period.The blue and orange dashed lines correspond to the equatorial and polar averaged approximations respectively.The Radial flow measured on the y axis is given by the dimensionless quantity dr/dt for each of the three flow equations.

Figure 3 .
Figure 3. Left Image: Plot of the maximum relative error between the equatorial and true flow over the range (r + , r I ).Right Image: Plot of the maximum relative error between the polar averaged and true flow over the range (r + , r I ).Both plots give the error over the entire parameter space (a, x inc ).

Figure 4 .
Figure 4. Orbital plots of plunging geodesics which asymptote to the ISSO in the infinite past with (a, r I ) = (0.9, 2.6) and θ = arccos(z).Here the cyan and orange planes are the azimuthally co-precessing and poloidally co-rotating planes respectively.Examples of the projection on these planes are shown in Fig. 5.

Figure 6 .
Figure 6.Orbital plot of a generic plunging geodesic with parameter values (a, E, L, Q) = (0.9, 0.94, 0.1, 12) and θ = arccos(z).The larger black sphere gives the horizon of the black hole where as the smaller sphere simply gives a point along the geodesics.Here the cyan and orange planes are the azimuthally co-precessing and poloidally co-rotating planes respectively.

Figure 7 .
Figure 7. Co-rotating orbital plots of a generic plunging geodesic with parameter values (a, E, L, Q) = (0.9, 0.94, 0.1, 12).Left Image: orthogonal projection of generic plunge onto the co-precessing azimuthal plane.Right Image: orthogonal projection of generic plunge onto co-rotating polar plane.The cyan and orange frames represent the projection onto the planes seen in 6 as they track a particle following the geodesic trajectory.