Finding high-order analytic post-Newtonian parameters from a high-precision numerical self-force calculation
aa r X i v : . [ g r- q c ] M a r Finding high-order analytic post-Newtonian parameters from a high-precision numerical self-force calculation
Abhay G. Shah,
1, 2, ∗ John L. Friedman, † and Bernard F. Whiting
4, 5, ‡ Dept of Particle Physics & Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel School of Mathematics, University of Southampton, Southampton SO17 1BJ, United Kingdom Center for Gravitation and Cosmology, Department of Physics,University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee, Wisconsin 53201, USA Department of Physics, P.O. Box 118440, University of Florida, Gainesville, Florida 32611-8440 G R ε C O , Institut d’Astrophysique de Paris — UMR 7095 du CNRS,Universit´e Pierre & Marie Curie, 98 bis boulevard Arago, 75014 Paris, France We present a novel analytic extraction of high-order post-Newtonian (pN) parameters that governquasi-circular binary systems. Coefficients in the pN expansion of the energy of a binary system canbe found from corresponding coefficients in an extreme-mass-ratio inspiral (EMRI) computation ofthe change ∆ U in the redshift factor of a circular orbit at fixed angular velocity. Remarkably, bycomputing this essentially gauge-invariant quantity to accuracy greater than one part in 10 , andby assuming that a subset of pN coefficients are rational numbers or products of π and a rational,we obtain the exact analytic coefficients. We find the previously unexpected result that the post-Newtonian expansion of ∆ U (and of the change ∆Ω in the angular velocity at fixed redshift factor)have conservative terms at half-integral pN order beginning with a 5.5 pN term. This implies theexistence of a corresponding 5.5 pN term in the expansion of the energy of a binary system.Coefficients in the pN series that do not belong to the subset just described are obtained toaccuracy better than 1 part in 10 − n at n th pN order. We work in a radiation gauge, findingthe radiative part of the metric perturbation from the gauge-invariant Weyl scalar ψ via a Hertzpotential. We use mode-sum renormalization, and find high-order renormalization coefficients bymatching a series in L = ℓ +1 / L behavior of the expression for ∆ U . The non-radiativeparts of the perturbed metric associated with changes in mass and angular momentum are calculatedin the Schwarzschild gauge. ∗ [email protected] † [email protected] ‡ [email protected]fl.edu I. INTRODUCTION
The principal approximation methods used to compute the inspiral of compact binary systems are the post-Newtonian expansion, in which an orbital angular velocity M Ω serves as the expansion parameter; and the self-forceor extreme-mass-ratio-inspiral (EMRI) approach, in which the small parameter is the mass ratio m /M of the binary’stwo components. Previous work by Blanchet et al. [1, 2] used an overlapping regime where both approximations arevalid to check the consistency of the renormalization methods used in the two approaches and to find numerical valuesof pN coefficients at orders beyond the reach of current analytical work.In the present paper, by working with much higher numerical accuracy – maintaining precision of at least one partin 10 in an EMRI computation of the perturbed orbital frequency and redshift factor, and by considering orbitsat much larger separation – with orbital radii extending to 10 M , we obtain two surprising results not seen in theprevious study: • (1) A subset of the pN parameters in lower-order analytical work had been found to be either rationals m/n orto be sums of rationals multiplied by powers of π , the Euler constant γ and square roots of integers. Our highprecision allows us to extract the exact analytical form of the subset of coefficients that are rationals or productsof the form rational × π from our numerical values up to 10 pN order, corresponding to corrections smaller by( v/c ) that the Newtonian value. • (2) In a pN expansion, conservative terms (terms even under the interchange of outgoing and ingoing radiation)are initially encountered at integral pN orders; dissipative terms (odd under the interchange of outgoing andingoing) first enter at 2.5 pN order. At higher order, dissipative terms can occur at either integral or half-integral order, depending on the details [3], while conservative terms enter at each integral order. We find thatconservative terms of 5.5 pN order appear in the expression for the redshift at fixed angular velocity (and thusin the expressions for the angular velocity at fixed redshift and in the expression for the energy of an orbit withgiven angular velocity). These quantities are conservative, and the presence of 5.5 pN terms was unexpected.The work reported here involves a binary system that, at zeroth order in the mass ratio, is described by a testparticle in a circular geodesic about a Schwarzschild black hole. At first order in m /M , the orbital parameters arealtered by a metric perturbation h αβ produced by the orbiting particle: The perturbed motion can be described bysaying that the particle moves on a circular geodesic of the metric g αβ + h ren αβ , where h ren αβ is the renormalized metricperturbation. The perturbed spacetime is helically symmetric, with a helical Killing vector k α that is tangent to theparticle’s 4-velocity, u α = U k α . (1)The constant of proportionality U is termed the redshift factor (first introduced by Detweiler [4]), and can be thoughtof as a contribution to the redshift, measured from the perturbed orbit of the mass m , that is independent of theinternal geometry of the mass. With the perturbed spacetime chosen so that the perturbed and unperturbed helicalKilling vectors coincide, the change in ∆ U at fixed angular velocity Ω has the form∆ U = − U h ren αβ u α u β =: − U H ren , (2)and it is invariant under gauge transformations generated by helically symmetric gauge vectors.A pN expansion of ∆ U , written in terms of a dimensionless R := ( M Ω) − / , has the form∆ U = − R + X n =1 α n R n +1 + X n =4 β n log RR n +1 + X n =7 γ n log RR n +1 + X n =10 δ n log RR n +1 + · · · , (3)where the post-Newtonian order n can take half-integral as well as integral values, starting at α . and β . . Thatis, integral values of coefficients of log k R/R n +1 start at pN order n = 3 k + 1; half-integral values appear to start at n = 3 k + 5 .
5, but we do not carry our numerical expansion far enough to find the first half-integral value for k = 2( γ . ) or for larger k . We compute ∆ U at a set of radii extending to 10 M and match to a series of this form. Asnoted in the abstract, the high numerical accuracy of ∆ U ( r ) allows us to find the coefficients α n , β n , and γ n witha precision at least as high as one part in 10 − n . At each pN order, we find that the coefficient of the highestoccurring power of log R is rational when n is an integer; and it has the form rational × π when n is a half-integer.The remaining coefficients for a given value of n are not of this form.Because the presence of α . and higher-order half-integral coefficients was not expected, we performed an elaborateset of checks. Our calculations were carried out in a radiation gauge, but, we repeated the entire numerical calculationof ∆ U in a Regge-Wheeler gauge, obtaining numerical agreement to 368 places of accuracy, that is, the retarded valuesof h uu for each ℓ -mode in the RG and RWZ gauges agree to more than 368 digits). This serves as a demanding testof both the numerical code and of the analytical computation on which it is based. Because the numerical calculationis performed in Mathematica, the comparison is also a check of Mathematica’s claimed numerical precision. AdrianOttewill and Marc Casals kindly used their codes to perform an independent radiation-gauge computation to comparewith ours at double-precision accuracy for small R. Specifically, for the s = ℓ = m = 2 term, we compared our valuesof the invariant, A lm R H R ∞ (see Eq. (7) below), at r/M = 10 , . Finally, we analytically computed α . (seeSect. III A).In Sect. II we briefly review the calculation of the renormalized ∆ U in a modified radiation gauge. In Sect. III wepresent the results of matching a sequence of values ∆ U ( r ) to a series of the form (3).We work in gravitational units ( G = c = 1) and use signature + − −− to conform to Newman-Penrose conventions. II. REVIEW OF ∆ U COMPUTATION
We consider a particle of mass m orbiting a Schwarzschild black hole of mass M . At zeroth order in m /M , thetrajectory is a circular orbit. In Schwarzschild coordinates, its angular velocity is Ω = p M/r , and its 4-velocity isgiven by u α = U ( t α + Ω φ α ) , with U = u t = 1 p − M/r . (4)We compute the change ∆ U at first-order in m /M in a modified radiation gauge, as detailed in [5]. We briefly reviewthe formalism here, noting first that Eq. (3) for ∆ U involves a single component H ren of the renormalized metricperturbation.For multipoles with ℓ ≥
2, the metric perturbation can be found in a radiation gauge from the the spin-2 retardedWeyl scalar, ψ , which has the form [5–7], ψ ( x ) = ψ (0)0 + ψ (1)0 + ψ (2)0 , (5)with ψ (0)0 = 4 π m u t ∆ r X ℓm A ℓm [( ℓ − ℓ ( ℓ + 1)( ℓ + 2)] / R H ( r < ) R ∞ ( r > ) Y ℓm ( θ, φ ) ¯ Y ℓm (cid:16) π , Ω t (cid:17) , (6a) ψ (1)0 = 8 πi m Ω u t ∆ X ℓm A ℓm [( ℓ − ℓ + 2)] / Y ℓm ( θ, φ ) ¯ Y ℓm (cid:16) π , Ω t (cid:17) × n [ im Ω r + 2 r ] R H ( r < ) R ∞ ( r > ) + ∆ [ R ′ H ( r ) R ∞ ( r ) θ ( r − r ) + R H ( r ) R ′∞ ( r ) θ ( r − r )] o , (6b) ψ (2)0 = − π m Ω u t X ℓm A ℓm Y ℓm ( θ, φ ) ¯ Y ℓm (cid:16) π , Ω t (cid:17) × (cid:26) [30 r − M r + 48 M r − m Ω r − − r ( r − M ) + 6 im Ω r ( r − M )] R H ( r < ) R ∞ ( r > )+2(6 r − M r + 16 M r − r ∆ + im Ω∆ r )[ R ′ H ( r ) R ∞ ( r ) θ ( r − r ) + R ′∞ ( r ) R H ( r ) θ ( r − r )] re + r ∆ [ R ′′ H ( r ) R ∞ ( r ) θ ( r − r ) + R ′′∞ ( r ) R H ( r ) θ ( r − r ) + W[ R H ( r ) , R ∞ ( r )] δ ( r − r )] ) , (6c)where ∆ = r − M r ; the functions R H and R ∞ (indices ℓ, m are suppressed) are the solutions to the homogenousradial Teukolsky equation that are ingoing and outgoing at the future event horizon and null infinity, respectively,and a prime denotes their derivative with respect to r ; W[ R H ( r ) , R ∞ ( r )] = R H R ′∞ − R ∞ R ′ H ; and the quantities A ℓm ,given by A ℓm := 1∆ W[ R H ( r ) , R ∞ ( r )] , (7)are constants, independent of r . The functions R H and R ∞ are calculated to more than 350 digits of accuracy usingexpansions in terms of hypergeometric functions given in [8], namely R H = e iǫx ( − x ) − − iǫ ∞ X n = −∞ a n F ( n + ν + 1 − iǫ, − n − ν − iǫ, − − iǫ ; x ) , (8) R ∞ = e iz z ν − ∞ X n = −∞ ( − z ) n b n U ( n + ν + 3 − iǫ, n + 2 ν + 2; − iz ) . (9)where x = 1 − r M , ǫ = 2 M m
Ω and z = − ǫx . We refer the reader to [8, 9] for the derivation of ν (the renormalizedangular momentum), and the coefficients a n and b n . Here F and U are the hypergeometric and the (Tricomi’s)confluent hypergeometric functions.The computation of the spin-weighted spherical harmonics s Y ℓ,m ( θ, φ ) is done analytically using [7].Once ψ is computed, the components of the metric perturbation are found from Hertz potential, Ψ, whose angularharmonics are related to those of ψ by an algebraic equation,Ψ ℓm = 8 ( − m ( ℓ + 2)( ℓ + 1) ℓ ( ℓ −
1) ¯ ψ ℓ, − m + 12 imM Ω ψ ℓm [( ℓ + 2)( ℓ + 1) ℓ ( ℓ − + 144 m M Ω (10)where Ψ = P ℓ,m Ψ ℓm ( r ) Y ℓm ( θ, φ ) e − im Ω t and ψ = P ℓ,m ψ ℓm ( r ) Y ℓm ( θ, φ ) e − im Ω t . The components of the metricalong the Kinnersley tetrad are h = r ð Ψ + ð Ψ) , (11) h = r (cid:20) ∂ t − f ∂ t ∂ r + f ∂ r − r − M )2 r ∂ t + f (3 r − M )2 r ∂ r + r − M r (cid:21) Ψ , (12) h = − r √ (cid:18) ∂ t − f ∂ r − r (cid:19) ¯ ð Ψ , (13)where f = ∆ /r and the operators ð and ¯ ð , acting on a spin-s quantity, η , are given by ð η = − ( ∂ θ + i csc θ∂ φ − s cot θ ) η, ¯ ð η = − ( ∂ θ − i csc θ∂ φ + s cot θ ) η. (14)The metric recovered from ψ ren0 above only specifies the radiative part of the perturbations ( ℓ ≥
2) and the fullmetric reconstruction requires one to take into account the change in mass and angular momentum of the Schwarzschildmetric and are associated with ℓ = 0 and ℓ = 1 harmonics, respectively. The contribution to the full H from thechange in mass ( H δM ) and angular momentum ( H δJ ) of the Schwarzschild metric are given by (see Eqs. (137, 138)of [5]) H δM = m ( r − M ) r / ( r − M ) / , (15) H δJ = − M m r / ( r − M ) / . (16)The renormalization of H is described in detail in [5–7]. The related quantity ∆Ω that gives the O ( m ) change inthe the angular velocity of a trajectory at fixed redshift factor is∆Ω = − u φ u t H ren = ∆ Uu φ u t . (17) III. RESULTS
In this section we present the pN-coefficients of ∆ U . Prior to this work, the following analytical coefficients wereknown [1, 10]:∆ U = − R + − R + − R + − π R + − − γ + 10155 π − R + 64 log( R )5 R + −
956 log( R )105 R (18)We calculate ∆ U for a set of R -values from 1 × , × , × , × to 10 in logarithmic intervals of 10 with an accuracyof one part in 10 for R = 10 , 10 for R = 10 and 10 for R = 10 . We then match this data to a pN-seriesto extract the unknown coefficients. In doing so, we find non-zero half-integer ( n.
5) pN coefficients that come fromthe tail-of-tail terms in pN-computations [11]. To confirm its presence we analytically calculated the 5.5pN term – thecoefficient of 1 /R . and found that it agreed with the numerically extracted coefficient to 113 significant digits. (Theanalytic calculation is described briefly below.) The high accuracy of the numerically extracted coefficients, however,allows us to extract their exact analytical expressions, without an analytic calculation . For example, the numericallyextracted value of the 6-pN log-term is − . · · · (19)More than five repetition cycles of the string 398589065255731922 tells us that it is the rational number − / U analytically known = − R + − R + − R + − π R + − − γ + 10155 π − R + 64 log( R )5 R + −
956 log( R )105 R + − π R . + − R )567 R + 81077 π R . + 27392 log ( R )525 R + 82561159 π R . + − ( R )2205 R + − π log( R )55125 R . + − ( R )9823275 R + 99186502 π log( R )1157625 R . + 23447552 log ( R )165375 R . (20)A rational number with fewer than ten digits in its numerator and in its denominator is determined by the first elevendigits in its decimal expansion; thus if one assumes that the rationals occurring in the coefficients of (20) have thischaracter, they are uniquely determined by the numerical accuracy. Without the assumption, the probability that thefirst n digits in a decimal representation of a randomly chosen number will match a rational with n n and n d digits innumerator and denominator is less than 10 n n + n d − n , when n > n n + n d .After using the above analytical coefficients, a numerical fit for the other numerical coefficients in (3) gives thevalues listed in Table I. A. Exact, Analytic 5.5pN value
As mentioned above, as a check on the work, we analytically compute the 5.5pN term. To do so, we use the fact thatthe renormalization parameters that characterize the singular part of H ret have no n.5pN terms: The pN expansionof H sing does not include half-integer powers of 1 /R . Studying the pN-expansion of the first few multipoles of H ret ,we find that the 5.5pN term comes only from the ℓ = 2 , m = ± H ret . That is, the numerical coefficientof the 5.5-pN term we obtain by matching H ren coincides exactly with that obtained by matching the sum of the ℓ = 2 , m = ± H ret to a pN-series. The analytic calculation was thus restricted to the ℓ = 2 , m = ± R ( p ) H R ( q ) ∞ / ( R H R ′∞ − R ′ H R ∞ ) (where p and q , the number of radial derivates, each run from 0 to 2). Weuse the hypergeometric series Eqs. (8) and (9) to express each of these functions as Taylor series in powers of 1 /R .From these series, we obtain in turn the pN-expansions of the l = 2 , m = ± ψ , Ψ and their firsttwo radial derivatives and, finally, the pN-series of H ret2 , ± . IV. NUMERICAL EXTRACTION AND ERROR ANALYSIS
We describe in this section the way we numerically extract the pN-coefficients and check the accuracy with whichthey are determined. We compute ∆ U ( R ) for R = 1 × , × , × , × to 10 in logarithmic intervals of 10. Fromthis data, after subtracting the known terms of Eq. (18), we match it to X n =5 α n R n +1 + X n =6 β n log RR n +1 + X n =7 γ n log RR n +1 + X n =10 δ n log RR n +1 + · · · , . (21) Coefficient Numerical value α -243.1768144646743075872935889669380023473727281723278653952886830882794813055787844008820951887564926056965827710452637773038028704808 a α -1305.001381078709655741090068271713685159580884739476033307892025133498776905927112179825227138960576902431854 a α -6343.8744531990306527270512066053061390446046295187692031581328657892063930482892366 α -11903.4729472013044159758685624140826902285745341620173222629 α . -8301.37370829085581136384718573193317705504946743 α -32239.6275950925564123677060345920962 α . -10864.625586706244075245767 α -221316.52514302 α . × β a β a β -3176.929181153969206392338832692666088 β -7358.271055677 β . γ a See
Note added at the end of the paper.
TABLE I. Numerical values of the coefficients in the expansion (2) of ∆ U for which analytic expressions could not be inferred. The accuracy with which the coefficients are extracted depends on the number of terms in the series. The fit is donein
M athematica and the number of terms in Eq. (21) chosen to maximize the accuracy of the extracted coefficientsis calculated as follows. Since, the coefficient extracted depends on the number of terms (say k ) in the fitting series,we give another index k to some n-pN term, say the first unknown coefficient 5.5-pN term, α k (we omit the pN-indexfor simplicity for now). We then look at the fractional difference of the | α k ± /α k − | vs k , and the k at whichthe fractional difference is minimum, we choose that coefficient. For further details we refer the reader to Sec (V)of [7] where a similar fitting is done. The fitting procedure is done twice here - first to extract the new analyticalpN-coefficients (the terms in Eq. (20) minus Eq. (18)), and then we subtract them from data to do another fit toextract the coefficients in Table I. V. DISCUSSION
In [12] it was established that a relation exists between, on the one hand, coefficients in the pN expansion of the red-shift variable and, on the other hand, coefficients in the expansion of the pN binding energy and angular momentumfor the binary system; for explicit results, see, for example, Eq. (2.50a-d) and (4.25a-d) in [12]. Subsequently, usingessentially Eqs. (2.40), (4.19) and (4.23) in [12], Le Tiec at al. (see [13]) transformed these relations to obtain, forarbitrary pN order, elegant expressions for the energy and angular momentum directly in terms of the self-force redshift variable and its first derivative; see, in particular, Eqs. (4a-b) in [13].Self-force extensions of the pN binding energy and angular momentum have been long sought after, since they wereknown to have the potential to contribute to the effective-one-body (EOB) formulation (see [14, 15]) of the binaryinspiral problem — mimicking, as far as possible, the reduced mass form of the Newtonian problem, but in a fullyfour dimensional, space-time setting. Thus, in a follow-up paper to [13], Barausse et al. [16] found a very compactresult, expressing the relevant EOB function directly in terms of the self force variable alone; see their Eq. (2.14) forthis important relation, subsequently also reported in [17].There is a very clear synergy between self-force results, and their applications in pN and EOB work, and knowledgeof our new results will have an immediate impact though the application of the relations discussed in the previoustwo paragraphs. Since the completion of our calculation, a corresponding computation has been performed to directlyevaluate the 5.5pN coefficient through conventional pN analysis, in which it is known to arise from a tails-of-tailscontribution. The ensuing result, as reported in a companion paper [11], is in exact agreement with the 5.5pN termin our Eq. (20).
Note : The works cited in this section express Ω as x / / ( M + m ) rather than our R / /M , and use z ( x ) = 1 /U ( R )as the red shift variable. The notation used throughout the rest of this paper was first introduced by Detweiler [4]. Note added:
At 1:15 pm (GMT) on December 9th, 2013 we received email notification from Nathan Kieran Johnson-McDaniel [18] that α could be represented as205680256 + 7342080 γ − π + 28968960 log(2) − . An equivalent result, and an exact expression for α , have subsequently appeared in [19]. It has since been possibleto show that our numerical results for β and β can be represented by β = 51637225195457375 − γ − β = 769841899153496621125 + 1080642205 γ + 536131211025 log(2) − . An explanation of these results and the methods used to obtain them will be discussed in a forthcoming paper [20].
ACKNOWLEDGMENTS
We are indebted to Alexandre Le Tiec for pointing out [10] for the value of 4-pN coefficient. This work wassupported by NSF Grant PHY 1001515 to UWM, PHY 0855503 to UF, European Research Council Starting GrantNo. 202996 to WIS and the European Research Council under the European Unions Seventh Framework Programme(FP7/2007-2013)/ERC grant agreement no. 304978 to UoS. BFW acknowledges sabbatical support from the CNRSthrough the IAP, where part of this work was carried out. [1] L. Blanchet, S. Detweiler, A. L. Tiec, and B. F. Whiting, Phys Rev D, , 084033 (2010).[2] L. Blanchet, S. Detweiler, A. Le Tiec, and B. F. Whiting, in Mass and Motion in General Relativity , edited by L. Blanchet,A. Spallicci, and B. Whiting (2011) pp. 415–442.[3] L. Blanchet and T. Damour, Phys. Rev. D, (1988).[4] S. Detweiler, Phys. Rev. D, , 124026 (2008), arXiv:0804.3529 [gr-qc].[5] A. G. Shah, J. L. Friedman, and T. S. Keidl, Phys Rev D, , 084059 (2012).[6] T. S. Keidl, A. G. Shah, J. L. Friedman, D.-H. Kim, and L. R. Price, Phys. Rev. D, , 124012 (2010).[7] A. G. Shah, T. S. Keidl, J. L. Friedman, D.-H. Kim, and L. R. Price, Phys. Rev. D, , 064018 (2011).[8] S. Mano, H. Suzuki, and E. Takasugi, Prog.Theor.Phys., , 1079 (1996).[9] M. Sasaki and H. Tagoshi, Living Rev. Relativity, (2003).[10] D. Bini and T. Damour, Phys Rev D, , 121501(R) (2013).[11] L. Blanchet, G. Faye, and B. F. Whiting, “Half-integral conservative post-newtonian approximations in the redshiftobservable of black hole binaries,” (2013), arXiv:1312.2975[gr-qc].[12] A. L. Tiec, L. Blanchet, and B. F. Whiting, Phys Rev D, , 064039 (2012).[13] A. Le Tiec, E. Barausse, and A. Buonanno, Phys.Rev.Lett., , 131103 (2012), arXiv:1111.5609 [gr-qc].[14] A. Buonanno and T. Damour, Phys. Rev. D, (1999).[15] T. Damour, P. Jaranowski, and G. Schaefer, Phys. Rev. D, (2000).[16] E. Barausse, A. Buonanno, and A. Le Tiec, Phys.Rev., D85 , 064010 (2012), arXiv:1111.5610 [gr-qc].[17] S. Akcay, L. Barack, T. Damour, and N. Sago, Phys.Rev.,