{'fraction_non_alphanumeric': 0.03625251172136638, 'fraction_numerical': 0.011428332217012727, 'mean_word_length': 4.5158162087277764, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical high-order post-Newtonian expansions for spinning extreme mass ratio binaries

We present an analytic computation of Detweiler's redshift invariant for a point mass in a circular orbit around a Kerr black hole, giving results up to 8.5 post-Newtonian order while making no assumptions on the magnitude of the spin of the black hole. Our calculation is based on the functional series method of Mano, Suzuki and Takasugi, and employs a rigorous mode-sum regularization prescription based on the Detweiler-Whiting singular-regular decomposition. The approximations used in our approach are minimal; we use the standard self-force expansion to linear order in the mass ratio, and the standard post-Newtonian expansion in the separation of the binary. A key advantage of this approach is that it produces expressions that include contributions at all orders in the spin of the Kerr black hole. While this work applies the method to the specific case of Detweiler's redshift invariant, it can be readily extended to other gauge invariant quantities and to higher post-Newtonian orders.


I. INTRODUCTION
Binary black hole systems have been identified as one of the primary sources of gravitational waves for current and planned generations of gravitational-wave detectors [1][2][3]. Accurate models for the waveforms produced by gravitational-wave sources are a crucial component in the data-analysis pipeline used to extraction of information from gravitational wave observations. In the context of compact-object binary systems, the production of waveform models typically relies on one of three fundamental methods: numerical relativity (NR) simulations, post-Newtonian (PN) approximations, or gravitational selfforce (GSF) calculations using black hole perturbation theory.
Gravitational self-force calculations -involving, for example, a solar mass black hole or neutron star of mass m in orbit around a massive black hole of mass M -are based on a perturbative expansion of Einstein's equations, with the mass ratio m/M as a small expansion parameter. These so-called extreme mass ratio inspiral (EMRI) systems are well approximated by an expansion to linear order in m/M . An alternative approach to the two-body problem -valid when the constituents are far apart -is the post-Newtonian approximation, which expands the Einstein equations in v 2 /c 2 , where v is a representative velocity and c is the speed of light. In the context of binary systems, the post-Newtonian expansion maps onto an expansion in 1/r, where r is the separation of the two objects.
The problem of gauge freedom in general relativity is a constant source of difficulty in extracting meaningful information from calculations. It is often difficult to know how much of the difference between two different results is due to merely a difference in choice of gauge. The computation of gauge-invariant quantities provides a robust solution to this problem. Independent of the choice of gauge, the computation of gauge-invariant quantities should agree among different methods; any discrepancy between results can be confidently associated with errors in the method. Even better, a gauge-invariant quantity computed within one approach can often be used to inform -and even derive results within -other approaches to the same problem.
The strategy of making comparisons based on gauge invariant quantities has proven to be a fruitful method for cross-pollination of results between NR, PN, and GSF [4][5][6][7][8][9][10][11][12][13][14]. This work continues a programme [10,11,[15][16][17] to compute the gauge-invariant quantities that enable these valuable cross-comparisons. Our key new development is the incorporation of spin effects, by applying methods similar to those of [15] to the case of analytic GSF calculations in Kerr spacetime, allowing us to add important spin-dependent terms into PN and effectiveone-body (EOB) [18] models. Given that most (if not all) astrophysical black holes are expected to be spinning [19,20], these additional terms are crucial for faithfully representing the type of systems we expect gravitational wave detectors to observe.
A parallel effort by Bini, Damour and Geralico [21] has recently been successful in computing spin-dependent contributions to GSF-PN results -that is, approximations to the spacetime of a binary system using simultaneous expansions in the mass ratio and the binary separation -using a small-spin approximation. Our work provides two key advantages over these results: 1. The only expansion our method relies only is the standard PN expansion in the inverse separation, y ∼ 1/r (in addition to the standard expansion in mass-ratio used in all self-force calculations). Importantly, this allows to provide for the first time expressions that are exact in the spin, a, of the larger black hole. The validity of our expressions in the high-spin regime is particularly important given observational evidence for black holes with near-extremal (a > 0.95M ) spins [19,20].
2. Our regularization procedure is based on an independently obtained expression for spheroidalharmonic mode-sum regularization parameters derived from the Detweiler-Whiting singular field.
In addition, our results provide a valuable independent check; indeed, they identify what we believe to be an error in some terms in the order y 9.5 PN coefficient for ∆U given by Ref. [21]. This conclusion has been confirmed by an independent check against numerically derived values [22]. This paper follows the conventions of Misner, Thorne and Wheeler [23]; a "mostly positive" metric signature, (−, +, +, +), is used for the spacetime metric, the connection coefficients are defined by Γ λ µν = 1 2 g λσ (g σµ,ν + g σν,µ − g µν,σ ), the Riemann tensor is R α λµν = Γ α λν,µ − Γ α λµ,ν + Γ α σµ Γ σ λν − Γ α σν Γ σ λµ , the Ricci tensor and scalar are R αβ = R µ αµβ and R = R α α , and the Einstein equations are G αβ = R αβ − 1 2 g αβ R = 8πT αβ . Standard geometrized units are used, with c = G = 1, but we include the explicit dependence on G and c in our post-Newtonian expansions in cases where they are convenient for post-Newtonian order counting. We use the spherical Boyer-Lindquist coordinates {t, r, θ, φ} for the background Kerr spacetime and write tensors in terms of these coordinate components.

II. REDSHIFT FOR CIRCULAR GEODESICS IN KERR SPACETIME
A. Circular, equatorial geodesic orbits in Kerr spacetime In this work we are interested in the case of a point mass on a circular equatorial geodesic in Kerr spacetime. In Boyer-Lindquist coordinates, the line element for a Kerr black hole of mass M and spin parameter a is given by (2.2) Circular equatorial orbits can be parametrized by the orbital frequency, Ω, which is related to the Boyer-Lindquist radius of the orbit, r 0 , by Adopting the convention that u ϕ (and hence the orbital angular momentum) is always positive, the orbital angular momentum and energy per unit mass for such orbits are given by Within this convention, prograde and retrograde orbits are distinguished by the sign of a: a > 0 for the former; a < 0 for the latter.

B. Redshift invariant
The key result of this work is the computation of a post-Newtonian expansion of Detweiler's redshift invariant [24] for circular equatorial orbits in Kerr spacetime. In this context, the redshift invariant is the constant of proportionality between the particle's 4-velocity, u α , and the helical Killing vector of the system, k α , i.e. (2.6) In the background spacetime, U is equal to the time component of the 4-velocity, u t . Provided k α is a Killing vector of both the background and perturbed spacetime (i.e. they both share the same helical symmetry) the contribution to U from the linear order metric contribution is simply given by This is invariant in the sense that it does not change under helically symmetric gauge transformations that also respect the helical symmetry of the worldline. When making comparisons between calculations using different gauges, care must be taken to ensure that these criterion are satisfied, or that the appropriate transformation of ∆U is accounted for [25].

A. Perturbations of Kerr spacetime
In the Kerr spacetime Teukolsky showed, using the Newman-Penrose formalism, that the dynamics of a metric perturbation are described by the evolution of tetrad components of the Weyl tensor [26,27]. This requires choosing a set of four null vectors e α i , i = 1, ...4, two real and two a complex conjugate pair. Using Boyer-Lindquist coordinates, a particular set of null vectors -the Kinnersly tetrad -is Using this tetrad, five of the spin-coefficients vanish, leaving us with The perturbed Weyl tensor C αβγδ , then has two components of interest which both independently contain all of the radiative information about perturbations of Kerr spacetime. Teukolsky's key insight was that the equations for ψ 0 and ψ 4 decouple, and each satisfy a separable PDE now known as the Teukolsky equation: respectively. In terms of tetrad components of the stress tensor, T ij = T µν e i µ e j ν , the source terms for the Teukolsky equation are where D = l µ ∂ µ , ∆ = n µ ∂ µ and δ = m µ ∂ µ . In this work we are interested in a perturbation sourced by a point particle, where a subscript-0 denotes evaluation on the worldline. In particular, for circular orbits θ 0 = π 2 , ϕ 0 = Ωt, so that T ij are the tetrad components of the stress tensor with the tetrad vectors evaluated at the position of the particle. The Teukolsky master equation, Eq. Here, s S ℓm (θ, ϕ; aω) are the oblate spin-weighted spheroidal harmonics with defining equation 1 sin θ d dθ sin θ d dθ − a 2 ω 2 sin 2 θ − (m + s cos θ) 2 sin 2 θ − 2aωs cos θ + s + 2maω + s λ ℓm s S ℓm (θ, ϕ; aω) = 0. (3.8) The spin-weighted spheroidal harmonics are orthonormal on the two-sphere, so that in the limit aω → 0 they coincide with the standard spin-weighted spherical harmonics s Y ℓm (θ, ϕ) (this normalisation is consistent with the Meixner-Schäfke convention for the spheroidal Legendre functions,

S
(1) mℓ (aω, cos θ), defined in [28]). The radial functions are solutions to where K = (r 2 + a 2 )ω − am, and where the source term on the right hand side is obtained from a mode decomposition of T 0 or T 4 . In the next section, we will develop analytic expressions for the appropriate homogeneous solutions to this equation.

B. Homogeneous Solutions of the Teukolsky equation
We now seek solutions to Eq. (3.10) for s = ±2, ℓ ≥ |s|, −ℓ ≤ m ≤ ℓ and ω ∈ R. To this end, we build up a complete set of ℓ modes from three distinct sections. For all ℓ ≥ 2 we can use the solutions to the Teukolsky equation give by Mano, Suzuki and Takasugi, as detailed in a review article by Sasaki and Tagoshi [29]. Furthermore, for sufficiently large values of ℓ the regularity of these solutions can be exploited with an ansatz to produce general expressions for arbitrary (large) ℓ and m. As such, in practice we: (i) compute a finite number of specific-ℓ values (with the exact number of required specific-ℓ values needed depending on the final PN order desired) using the MST solutions; and (ii) obtain the remaining modes from our ansatz. The third section of the solution relates to the non-radiative modes, and is addressed in Sec. III D.

Low ℓ modes: MST
Descriptions of the construction of the MST series solutions to the Teukolsky equation are widely available in the literature [29], so we will give here only a brief overview.
In computing PN expansions, we find two natural small parameters; the frequency ω = mΩ and the inverse of the radius, 1/r, which is related to the orbital frequency, Ω, via Eq. (2.3). This double expansion is handled indirectly by expanding in the inverse of the speed of light, 1/c, and introducing the auxiliary variables X 1 = GM/r and X 1/2 2 = ωr. Dimensional analysis shows that each of these variables carry a factor of 1/c, and are of the same order in the large-r limit. Towards the end of our calculations we will change back to a single expansion variable, namely y = (M Ω) 2/3 and use this relation to define X 1 (y) and X 2 (y).
The homogeneous solution satisfying retarded boundary conditions at the horizon can be written as a convergent infinite sum of hypergeometric functions Here ν is the well-known "renormalised angular momentum"; for the purpose of this section the critical feature is that and more convenient expressions for they hypergeometric function can be obtained by using the decomposition We denote the first and second terms on the right hand side by F 1 and F 2 , respectively; this leads to the R in = High order expansions of these expressions can prove computationally quite expensive if done without care. An examination of the leading order behaviour in η = 1/c of each term in the sums involved in both F 1 and F 2 helps to minimise the expansion of each n-value and determine a look-up table for truncating the infinite sum to a given order in η. In doing this one must take care of the irregular behaviour of the the series coefficients a ν n , the Γ functions, and the 2 F 1 with changing n. We summarise our results in Tables I, II and III. One minor complication we find is that the static m = 0 case must be handled separately. In the circular orbit case this problem can be circumvented; when m = 0 then also ω = mΩ = 0. Closed form analytic expressions for the solutions in this case were derived by Teukolsky [30], and more recently given in [31] and [32].
The homogeneous solution satisfying radiative boundary conditions at infinity can be written as a sum over irregular confluent hypergeometric functions with the same series coefficients Then, Note that the prefactors here have been taken to agree with Sasaki and Tagoshi [29]. Now using the standard identity we may split (3.13) into two more manageable pieces (3.14)Ũ 2 (n + 1 + s + ν − iǫ, 2n + 2 + 2ν; −2iz) = Γ (2n + 1 + 2ν) The leading order behaviour in η of each term in the sums involved in bothŨ 1 andŨ 2 are summarised in Tables IV and V. Once again we treat the m = 0 case differently. The tables are presented, however the analytic expressions given in [32] can alternatively be used.

Phase extraction and the large ℓ modes: A PN ansatz
As described in [33], for sufficiently large values of ℓ, in the Schwarzschild case the MST solutions to the Regge-Wheeler equation, X in/up can be written in the form: where ψ in and ψ up are r-independent phase factors and the A i 's and B i 's are pure polynomials in X 1 , X 2 1/2 . (For low-ℓ values, this expression is corrupted by logarithms and odd powers of η.) At this point, two further optimisations help to dramatically improve the efficiency of calculations. In constructing the retarded Green function, the phase factors will drop out since they amount to an irrelevant normalisation. Throwing them away and working with the homogeneous solutions without the phase factors leaves expressions which are orders of magnitude smaller in complexity. A further simplification can be found by starting from the a large ℓ ansatz and solving for the a (i,j) by demanding we have a solution of the Regge-Wheeler equation. This provides an efficient way to generate homogeneous solutions with ℓ, m left unspecified, avoiding the complexity of the MST solutions.
We find that the situation with the Teukolsky equation 3(|n| + 1) 3ℓ + 6 3ℓ + 3 3|n| The leading behaviour of the hypergeometric functions appearing as the first and second terms in Eq. (3.12). The terms marked with a † highlight the interplay between parameters and argument; see [33] for a full discussion.
is similar. Here the solutions can be re-expressed as where once again up until an ℓ-dependent power of η the A i and B i are pure polynomials in X 1 , X 2 1/2 and the essentially irrelevant phase functions remove significant complexity. As in the Regge-Wheeler case, we proceed by using this as an ansatz for solutions with ℓ, m unspecified and obtain general expressions for the large-ℓ homogeneous solutions.

Teukolsky-Starobinsky identities
In the previous section, the choice of spin, s, was arbitrary; s = +2 corresponds to ψ 0 while s = −2 corresponds to ψ 4 . Only one or the other is required in order to obtain the full radiative metric perturbation. However, in some situations it may be more convenient to have one or the other. Fortunately, rather than repeating a lengthy calculation for both there is a convenient shortcut. Given our set of spin −2 homogeneous solutions, we can easily calculate the spin +2 solution via a set of differential transformations known as the Teukolsky-Starobinsky identities. These are discussed in detail in many places, for example see Ref. [34]; for completeness we repeat the final result here. Given either an "in" or It is worth noting that on the computational side, upon doing this transformation, we will reintroduce rindependent terms of the form X 1 X 2 1/2 that can once again be extracted as components in the phase and essentially ignored.

C. Reconstructed metric perturbation in radiation gauge
The procedure for building up the components of the metric perturbation from the Weyl scalars [35][36][37][38][39] involves the construction of a Hertz potential, Ψ, from which gives the metric perturbation can be computed using where the overall minus sign here accounts for our mostly positive sign convention for the metric. For simplicity, from this point forward we restricting ourselves to the outgoing radiation gauge (ORG) defined by n µ h µν = 0.
In this gauge the Hertz potential itself is a solution of the spin-2 Teukolsky equation. Calculating Ψ can be done either using ψ 0 or ψ 4 in a variety of ways. For example, in [32] expressions are given for the Ψ in terms of the asymptotic amplitudes of ψ 4 , whereas in [40] it is constructed by inverting a differential operator which simplifies in the circular orbit case. Therefore to proceed one must solve either the s = 2 or s = −2 Teukolsky equation. However from a practical point of view it is not particularly important which is chosen, as the Teukolsky-Starobinsky identities can be used to transform between the two.
In this work we chose to construct ψ 0 using the s = 2 homogeneous solutions. To do this we construct the retarded Green function

D. Metric completion
It is well known that a metric reconstruction procedure based on ψ 0 or ψ 4 does not yield the whole radiationgauge metric perturbation. Wald [41] showed that for a Kerr background the remaining part of the metric perturbation can be fully attributed to perturbations to the mass and angular momentum of the background black hole (often, these are informally described as "ℓ = 0 and ℓ = 1 parts" of the perturbation). Here, we follow the standard procedure [42] and incorporate this contribution using analytic expressions. In particular, defining H = 1 2 h αβ u α u β , the additional contributions corresponding to perturbations to the mass and angular momentum are H δM = There is a subtlety here in that these contributions are not smooth on the worldline, and this non-smoothness could introduce additional contributions to the regularization procedure (see [43] for a more detailed discussion). This is a complicated issue worthy of a detailed independent analysis; here, we merely follow the standard procedure of evaluating the contributions from the completion part in the limit r → r + 0 and use the mode-sum regularization procedure described in the next section.

E. Regularization
We adopt a variation of the standard mode-sum regularization approach in order to extract a finite value from the divergent retarded metric perturbation. Traditionally, this mode-sum approach is written in terms of a sum over regularized spherical harmonic modes. This has the distinct disadvantage of requiring a cumbersome projection of the modes of the retarded metric perturbation onto scalar spherical harmonics. We have avoided this unnecessary step by instead deriving a mode-sum formula for the spheroidal harmonic modes that naturally arise from solutions of the Teukolsky equation. The full details of this derivation will be given in a forthcoming work; here we merely highlight the key results.
The derivation of our mode-sum formula is conceptually similar to previous derivations in terms of spherical harmonics, i.e.
1. Work in a spherical coordinate system, (α, β) in which the particle is instantaneously located at the north pole, α = 0.
2. Obtain a local coordinate expansion of the contribution to H from the Detweiler-Whiting singular field.
3. Decompose this coordinate approximation into spin-0 spheroidal harmonics, where the spheroidal harmonics are defined with respect to a coordinate system (θ, ϕ) in which the worldline is in the equatorial plane, θ = 0. The decomposition process makes use of the relation between the modes in the (α, β) coordinate system (where the mode decomposition can be most-easily done analytically) and the modes in the (θ, ϕ) coordinate system (where the retarded-field modes are most easily obtained).
4. Sum over m (azimuthal) modes to obtain a modesum formula.
The result of this process is a mode-sum formula for computing the regularized redshift invariant, where and with K(k) and E(k) being complete elliptic integrals. For the purposes of this work, we require the post-Newtonian expansion of this regularization parameter, which is given by

V. DISCUSSION
In this work, we have presented results for the linearin-mass ratio contribution to Detweiler's redshift invari-ant, ∆U , for a quasi-circular binary black hole system, in the case where the larger black hole is spinning. Our results are given as a PN-type expansion in the inverse separation, y, of the binary, but are otherwise exact. In particular, they are valid for astrophysically-important cases where the spin of the larger black hole is arbitrarily large.
There are several clear future directions for this work, in particular: • An extension to eccentric and inclined orbits would allow for a more complete exploration of the parameter space. The components for such a calculation are readily available -eccentric orbits have been studied without spin in [17, 45,46], while spin contributions have been considered here and in [21]so such a calculation would merely require an appropriate combination of the two approaches.
• An extension to higher PN orders. As the method we present here is totally algorithmic and implemented as a Mathematica code, it is straightforward (but more computationally expensive) to apply it to higher orders if the demand arises. Indeed, an online repository of our results [44] will be updated as future results become available.
• The application of our method to the calculation of other gauge invariant quantities, such as the spinprecession [10], tidal [11], and octupolar [16] invariants. These are obtained from the same metric perturbations used for the redshift invariant, so their calculation would be a straightforward application of the results presented here.
In addition to these, it is likely that an application of our results to improving and informing PN and EOB theories -for example by computing GSF contributions to the potentials appearing in EOB theory -would yield valuable improvements to both PN and EOB models.

ACKNOWLEDGMENTS
We thank Abhay Shah for performing numerical checks of our results. We also thank Abhay Shah, Bernard Whiting, Seth Hopper, Niels Warburton, Marc Casals and Sarp Akcay for many helpful discussions. B.W. was supported by the Irish Research Council, which is funded under the National Development Plan for Ireland. C.K. was funded under the Programme for Research in Third Level Institutions (PRTLI) Cycle 5 and co-funded under the European Regional Development Fund. This material is based upon work supported by the National Science Foundation under Grant Number 1417132.