Factorization and resummation: A new paradigm to improve gravitational wave amplitudes. III: the spinning test-body terms

TL;DR

This paper introduces new methods for resumming gravitational wave amplitudes from spinning test bodies around black holes, significantly improving accuracy over previous models by combining analytical and numerical techniques.

Contribution

It develops novel resummation procedures for gravitational wave amplitudes that outperform existing methods, based on detailed analytical and numerical calculations of fluxes and waveforms.

Findings

New resummation procedures improve waveform accuracy.

Analytical fluxes match numerical results within 10^{-5}.

Enhanced models suggest modifications for current waveform approaches.

Abstract

We present new calculations of the energy flux of a spinning test-body on circular orbits around a Schwarzschild black hole at linear order in the particle spin. We compute the multipolar fluxes up to $\ell=m=6$ using two independent numerical solvers of theTeukolsky equation, one in the time domain and the other in the frequency domain. After linearization in the spin of the particle, we obtain an excellent agreement ($\sim 10^{-5}$) between the two numerical results.The calculation of the multipolar fluxes is also performed analytically (up to $\ell=7$) using the post-Newtonian (PN) expansion of the Teukolsky equation solution; each mode is obtained at 5.5PN order beyond the corresponding leading-order contribution. From the analytical fluxes we obtain the PN-expanded analytical waveform amplitudes. These quantities are then resummed using new procedures either based on the factorization of the orbital contribution (and resumming it independently from the spin-dependent factor) or on the factorization of the tail contribution solely for odd-parity multipoles. We compare these prescriptions and the resummation procedure proposed in Pan et al. [Phys. Rev. D 83 (2011) 064003] to the numerical data. We find that the new procedures significantly improve over the existing one that, notably, is inconsistent with the numerical data for $\ell+m=\text{odd}$ multipoles already at low orbital frequencies. Our study suggests that the approach to waveform resummation used in current effective-one-body-based waveform models should be modified to improve its robustness and accuracy all over the binary parameter space.

Tables

Table 1. Table 1: Fractional differences between the resummed and the numerical h ^ ℓ ​ m subscript ^ ℎ ℓ 𝑚 \hat{h}_{\ell m} ’s at the LSO shown in Fig.  5 . The fractional difference is defined as | h ^ num − h ^ analyt | / h ^ num superscript ^ ℎ num superscript ^ ℎ analyt superscript ^ ℎ num |\hat{h}^{\rm num}-\hat{h}^{\rm analyt}|/\hat{h}^{\rm num} computed at x LSO subscript 𝑥 LSO x_{\rm LSO} .

$(\ell ,m)$	${\Delta \left(x\right)\hat{h}\left(x\right)/\hat{h}|}_{x=x_{\mathrm{LSO}}}$ for $-0.9\le \sigma \le +0.9$
	$-0.9$	$-0.5$	$+0.5$	$+0.9$
$(2,2)$	$8\times {10}^{-4}$	$5\times {10}^{-5}$	0.0023	0.0024
$(2,1)$	$4\times {10}^{-5}$	$7\times {10}^{-5}$	$5\times {10}^{-5}$	$4\times {10}^{-4}$
$(3,3)$	0.0035	0.0012	$2\times {10}^{-4}$	0.0056
$(3,2)$	$9\times {10}^{-4}$	$5\times {10}^{-4}$	0.0041	0.0086
$(3,1)$	0.0057	0.0018	0.0023	0.0091
$(4,4)$	0.0098	0.0046	0.0017	0.0085
$(4,3)$	0.0049	0.0029	0.013	0.029
$(4,2)$	0.0047	0.0013	0.0013	0.007
$(4,1)$	$7\times {10}^{-5}$	$3\times {10}^{-5}$	$4\times {10}^{-4}$	0.001
$(5,5)$	0.014	0.0055	0.0072	0.042
$(5,4)$	0.012	0.0068	0.023	0.05
$(5,3)$	0.0029	$6\times {10}^{-4}$	${10}^{-4}$	0.0036
$(5,2)$	$8\times {10}^{-4}$	$5\times {10}^{-4}$	0.018	0.046
$(5,1)$	0.019	0.0067	0.0097	0.032
$(6,6)$	0.02	0.0087	0.013	0.078

Equations

$$\frac{\mathrm{D}{p}^{\mu }}{d\tau }=$$

$$\frac{\mathrm{D}{S}^{\mu \nu }}{d\tau }=$$

$$\sigma \equiv \frac{s}{\mu M}.$$

$$\Omega$$

$$\mathcal{O}{}_{-2}\psi (t,r,\theta ,\varphi )=4\pi T,$$

$${}_{-2}\psi =\frac{1}{2\pi }\int e^{-i\omega t}{\psi }_{\ell m}(r){}_{-2}{S}_{\ell m}(\theta ,\varphi )d\omega ,$$

$$\mathcal{D}{\psi }_{\ell m}(r)=4\pi T_{\ell m}(r).$$

$$T_{\ell m}=$$

$$+A_3{\delta }^{\prime \prime \prime }(r-r_p).$$

$${\psi }_{\ell m}(r)=C_{\ell m}^{+}{R}_{\ell m}^{+}(r)+C_{\ell m}^{-}{R}_{\ell m}^{-}(r),$$

$$C_{\ell m}^{\pm }$$

$$=B_0{R}_{\ell m}^{\mp }(r_p)+B_1{R}_{\ell m}^{\mp }{}^{\prime }(r_p)+B_2{R}_{\ell m}^{\mp }{}^{\prime \prime }(r_p)$$

$$+B_3{R}_{\ell m}^{\mp }{}^{\prime \prime \prime }(r_p).$$

$$\frac{d{E}^{\infty }}{dt}=\sum _{\ell m}{F}_{\ell m}=\sum _{\ell m}\frac{|C_{\ell m}^{+}{|}^2}{4\pi m^2{\Omega }^2}$$

$$F_{\ell m}=F_{\ell m}^{(N,ϵ)}{\hat{F}}_{\ell m},$$

$$h_{\ell m}(x)=\mathcal{R}{h}_{\ell m}^{(N,ϵ)}(x){\hat{h}}_{\ell m}^{(ϵ)}(x),$$

$$x^{\mathrm{LSO}}=\frac{1}{6}+\frac{\sigma }{12\sqrt{6}}.$$

$$|{\hat{h}}_{\ell m}^{(ϵ)}|={\hat{S}}^{(ϵ)}|{\hat{h}}_{\ell m}^{\mathrm{tail}}(x)|({\rho }_{\ell m}{)}^{\ell },$$

$$|{\hat{h}}_{\ell m}^{\mathrm{tail}}(x){|}^2=\frac{4\pi m{x}^{3/2}{\prod }_{s=1}^{\ell }{\left(s^2+2m{x}^{3/2}\right)}^2}{(\ell !{)}^2\left(1-e^{-4\pi m}\right)}.$$

$$\hat{E}(x)$$

$$=\frac{1-2x}{\sqrt{1-3x}}-\frac{x^{5/2}\sigma }{\sqrt{1-3x}},$$

$$\hat{j}(x)$$

$$=\frac{1}{\sqrt{1-3x}}-\sqrt{x}\left(1-\frac{1-4x}{\sqrt{1-3x}}\right)\sigma ,$$

$$\rho_{\ell m}=T_{n}\Biggl{[}\frac{\sqrt{\hat{F}_{\ell m}}}{|\hat{h}^{\rm tail}_{\ell m}(x)|\hat{S}^{(\epsilon)}}\Biggr{]}^{1/\ell},$$

$${\rho }_{\ell m}={\rho }_{\ell m}^{\mathrm{orb}}+{\rho }_{\ell m}^{\sigma }.$$

$${\rho }_{\ell m}(x,\sigma )=P_d^n[{\rho }_{\ell m}^{\mathrm{orb}}]\overline{{\hat{\rho }}_{\ell m}^{\sigma }},$$

$$|{\hat{h}}_{\ell m}^{(0)}(x,\sigma )|=\hat{E}|{\hat{h}}_{\ell m}^{\mathrm{tail}}|{\left[{\rho }_{\ell m}(x,\sigma )\right]}^{\ell }.$$

$$|{\hat{h}}_{\ell m}^{(1)}(x,\sigma )|\equiv |{\hat{h}}_{\ell m}^{\mathrm{orb},(1)}|+\sigma |{\hat{h}}_{\ell m}^{\sigma ,(1)}|,$$

$${\rho }_{\ell m}^{\mathrm{orb}}=T_n\left[{\left(\frac{|{\hat{h}}_{\ell m}^{\mathrm{orb},(1)}|}{|{\hat{h}}_{\ell m}^{\mathrm{tail}}|{\hat{j}}^{\mathrm{orb}}}\right)}^{1/\ell }\right],$$

$$|\tilde{h}{|}_{\ell m}^{\sigma }=T_n\left[\frac{|{\hat{h}}_{\ell m}^{\sigma ,(1)}|}{|{\hat{h}}_{\ell m}^{\mathrm{tail}}|}\right].$$

$$|{\hat{h}}_{\ell m}^{(1)}(x,\sigma )|=|{\hat{h}}_{\ell m}^{\mathrm{tail}}|\left\{{\left({\rho }_{\ell m}^{\mathrm{orb}}\right)}^{\ell }{\hat{j}}^{\mathrm{orb}}+\sigma |{\tilde{h}}_{\ell m}^{\sigma ,(1)}|\right\}.$$

$$|{\hat{h}}_{\ell m}|={\hat{j}}^{\mathrm{orb}}|{\hat{h}}_{\ell m}^{\mathrm{tail}}|(P_n^m[{\rho }_{\ell m}^{\mathrm{orb}}]{)}^{\ell }(1+{\hat{f}}_{\ell m}^{\sigma }).$$

$${\hat{f}}_{21}^{\sigma }$$

$${\hat{f}}_{33}^{\sigma }$$

$${\hat{f}}_{31}^{\sigma }$$

Full text

I Introduction

Test-mass results have been crucial to devise robust resummation techniques for the truncated post-Newtonian expansion that give access to analytical gravitational waveform and fluxes for circularized, nonprecessing, binaries Damour et al. (1998, 2009); Pan et al. (2011a); Nagar and Shah (2016); Messina et al. (2018). Such resummed waveform, and related fluxes, are one of the crucial building blocks of effective-one-body (EOB) waveform models for coalescing relativistic binaries Nagar et al. (2019); Akcay et al. (2019); Nagar et al. (2018); Cotesta et al. (2018); Hinderer and Babak (2017); Bohé et al. (2017). Up to now, resummation of PN-expanded analytical result is the only approach that can be adopted to improve the behavior of the PN-expansions in the strong-field, fast velocity regime up to merger Damour and Nagar (2009). From the very beginning of this endeavor Damour et al. (1998) the development (and testing) of resummation techniques has been driven by comparisons between some analytically resummed waveform and numerical waveforms (or fluxes) generated by a nonspinning particle inspiralling and plunging into a Schwarzschild or a Kerr black hole Damour and Nagar (2007); Harms et al. (2014). By contrast, none of the resummation approaches routinely used in state-of-the-art EOB models Nagar et al. (2018); Bohé et al. (2017) has been tested in the special case where the particle (which models a test black-hole) is spinning. This was not done up to now for at least two reasons: (i) on the one hand, robust and accurate numerical computations of the energy fluxes from a spinning particle on circular orbits around a Kerr black became available only recently Lukes-Gerakopoulos et al. (2017); Harms et al. (2016a, b); (ii) on the other hand, the analytical PN knowledge of the fluxes of a spinning particle around a Kerr black hole was only known at global 2.5PN order Tanaka et al. (1996) and only recently pushed to 3.5PN accuracy Cotesta et al. (2018). This paper builds on previous works and improves them along two directions: (i) the numerical fluxes of Refs. Harms et al. (2016b, a) are recomputed, in the time-domain, at an improved accuracy and increasing the number of multipoles. In addition, they are compared with an analogous calculation performed with a completely independent numerical code in the frequency domain, finding excellent consistency between the two methods once the results are linearized in the particle spin; (ii) though we here only consider the case of a spinning particle around a Schwarzschild black hole, the 2.5PN accurate results of Ref. Tanaka et al. (1996) are pushed to much higher PN order, namely relative 5.5PN accuracy for all multipoles of the flux up to $\ell=7$. The availability of such new PN information, at high order, allows us to extensively test the standard waveform resummation techniques of Refs. Damour et al. (1998, 2009); Pan et al. (2011a) in a corner of the binary parameter space that had not been covered before. Similarly, we use these new numerical data to check new resummation approaches proposed recently in Refs. Nagar and Shah (2016); Messina et al. (2018) and that are going to be partly incorporated in the next generation of EOB waveform models Riemenschneider et al. .

The paper is organized as follows. In Sec. II we summarize the analytical and numerical approaches to compute the gravitational wave fluxes from a spinning particle around a Schwarzschild black hole. From the energy fluxes, decomposed in multipoles, we obtain the gravitational waveform amplitudes, both numerically and analytically. The aim of Sec. III is to compare the numerical waveform amplitudes with several analytical representations, either in PN-expanded form or using some resummation technique. Concluding remarks are reported in Sec. IV. The paper is completed by technical Appendixes that explicitly report the outcome of the PN calculations. Throughout this work we use geometrized units such that $G=c=1$. We also define $M$ as the mass of the primary black-hole and $\mu$ as the mass of the secondary black hole such that $M\gg\mu$.

II Energy fluxes emitted by a spinning particle around a Schwarzschild black hole

In this section we consider the energy flux radiated in gravitational waves by a spinning particle on a circular orbit of radius $r_{0}$ around a Schwarzschild black hole. We will restrict our attention to the case where the particle’s spin axis is aligned with the orbital angular momentum. We compute the radiated flux both analytically, via high-order post-Newtonian calculations, and numerically, using two independent approaches.

II.1 Post-Newtonian results

II.1.1 Dynamics of a spinning particle

The equations of motion of a spinning test body moving on a curved spacetime background (Schwarzschild) are given by the Mathisson-Papapetrou-Dixon equations (MPD) Mathisson (1937); Papapetrou (1951); Dixon (1970)

[TABLE]

where $p^{\mu}$ is the four-momentum, $u^{\mu}$ is the four-velocity ($\tau$ is the proper time), $R^{\mu}{}_{\nu\alpha\beta}$ is the Riemann tensor of the spacetime and $S^{\alpha\beta}$ is the spin-tensor. From the spin tensor we can define the spin magnitude, $s=(\frac{1}{2}S_{\mu\nu}S^{\mu\nu})^{1/2}$, from which we define the spin variable

[TABLE]

As in previous work Lukes-Gerakopoulos et al. (2017); Harms et al. (2016a, b), we consider $-1\leq\sigma\leq 1$. We will comment in Appendix A on the meaning of these limits and on the interpretation of the test-mass results presented here as part of the general case where the masses of the two bodies are comparable. The MPD equations do not form a closed systems of evolution equations and a closure, called spin supplementary condition (SSC), is required. The choice of a SSC amounts to choosing a centre of mass for the spinning body Costa and Natário (2015). For astrophysically relevant values of the spin, there are strong indications that any physically meaningful quantity will not depend upon the SSC choice Harms et al. (2016a); Lukes-Gerakopoulos et al. (2017). In this paper we choose the Tulczyjew-Dixon condition $S^{\mu\nu}p{}_{\nu}=0$ Tulczyjew (1965).

Throughout this section, we will work to linear order in $\sigma$ for which we have $p^{\alpha}=\mu u^{\alpha}+\mathcal{O}(\sigma^{2})$. The orbital frequency is then given by Tanaka et al. (1996)

[TABLE]

where $u=M/r_{0}$ and $r_{0}$ is the radius of a circular orbit expressed in Schwarzschild coordinates.

II.1.2 Calculation of the energy fluxes

In the Teukolsky approach, the problem of computing the energy and angular momentum fluxes at infinity from perturbations to either a Kerr or Schwarzschild black hole has been well laid out in the literature. Since the main aim of this work is to use the output of such calculations, we will only give an overview of the main ideas and refer the reader to Misao Sasaki and Hideyuki Tagoshi’s living review Sasaki and Tagoshi (2003) for an in-depth discussion of the topic.

We begin with the spin-weight $s=-2$ Teukolsky equation

[TABLE]

where $\mathcal{O}$ is a second-order differential operator and $T$ is formed from the stress energy of the perturbation $T_{\mu\nu}$. When working in the frequency domain, the problem is simplified greatly, and by separating ${}_{-2}\psi$ into radial and angular components, namely

[TABLE]

where ${}_{-2}S_{\ell m}(\theta,\phi)$ are the spin-weighted spheroidal harmonics which reduce to the spin-weighted spherical harmonics in the Schwarzschild limit. The radial functions $\psi_{\ell m}(r)$ then satisfy a second order ordinary differential equation

[TABLE]

In the case of a spinning particle on a circular orbit around a Schwarzschild black hole, the source term here takes the schematic form

[TABLE]

If we have a pair of homogeneous solutions satisfying retarded boundary conditions at the horizon and at infinity, $R_{\ell m}^{-}$ and $R_{\ell m}^{+}$ respectively, then the general solution to this equation is given by

[TABLE]

where

[TABLE]

The $B_{i}=B_{i}(r_{p},\sigma)$ are determined by completing the integral, $W$ is the invariant Wronskian of the two solutions and $\Delta=r(r-2M)$. Normalizing the homogeneous solutions so that $R_{+}\sim r^{3}e^{i\omega r^{*}}$ as $r\rightarrow\infty$, the energy flux emitted at radial infinity is given simply by

[TABLE]

where $\Omega$ is the orbital frequency given in Eq. (4).

We now briefly describe the calculation of the $F_{{\ell m}}$ functions appearing in Eq. (11) as a PN expansion, i.e an expansion for large orbital radius and small frequency. This calculation was first presented for a spinning particle on a circular orbit by Tanaka et al. Tanaka et al. (1996), though truncated at 2.5PN order. We will extend their calculation to a higher order in the expansion. To evaluate Eq. (10) we first need PN expansions for the homogeneous radial functions $R^{\pm}_{\ell m}$. The well established method for systematically computing these uses the solutions of Mano, Suzuki and Takasugi (MST)  Mano et al. (1996a, b). The MST solutions to the homogeneous radial Teukolsky equation are given as infinite series of either hypergeometric ${}_{2}F_{1}$ functions or irregular confluent hypergeometric $U$ functions. In the PN limit, restricting the frequencies to the allowed harmonics of the orbital frequency, these infinite series truncate at finite orders and, with modern algebraic software, are methodically Taylor expanded to a desired PN order. For in depth discussion of our calculations see Kavanagh et al. (2016). With these in hand the evaluation of Eq. (10) is straightforwardly accomplished by also Taylor expanding the $B_{i}$’s, thus giving PN expansions for the energy flux modes. At each order in the PN expansion we also work to linear order in the spin $\sigma$. In practice, we found it easy to consider modes up to $\ell=7$ and obtain each multipolar flux at 5.5PN accuracy beyond the leading-order contribution. More precisely, each PN-expanded multipolar contribution is factorized as

[TABLE]

where $F_{\ell m}^{(N,\epsilon)}$ is the Newtonian (or leading-order) prefactor, while $\hat{F}_{\ell m}$ is the PN correction. Here, $\epsilon=0,1$ depending on the parity of $\ell+m$. For each $(\ell,m)$, $\hat{F}_{\ell m}$ is given by a 5.5PN-accurate polynomial, i.e. it has the structure $1+x+x^{3/2}+x^{2}+x^{5/2}+\dots+x^{11/2}$, where $x=(GM\Omega/c^{3})^{2/3}={\cal O}(c^{-2})$ is the PN-ordering frequency parameter. The defining formulas for the $F_{\ell m}^{(N,\epsilon)}$ are given explicitly in Sec. B.1; the Newton-normalized PN-expanded multipolar fluxes, $\hat{F}_{\ell m}$, in Sec. B.2.

II.2 Numerical results

For computing numerically the radiated fluxes we employ two codes that solve the Teukolsky equation. One is a 2+1 time-domain code of Ref. Harms et al. (2016b, a) and the other a frequency domain code Akcay et al. . The details of each of these codes are presented elsewhere and so here we will only give a brief overview of each code and show that their results are consistent with one another.

II.2.1 Frequency domain approach

The method employed by the frequency domain numerical code follows closely the description of the PN calculation given in Sec. II.1.2 with the exception that the homogeneous solutions are computed numerically. This is done using the semi-analytic MST method (see Sasaki and Tagoshi (2003) for a review and Casals and Ottewill (2015); Casals et al. (2016) for extensions we use). A similar version of the code to compute the homogeneous solutions is publicly available as part of the Black Hole Perturbation Toolkit BHP . With the homogeneous solutions in hand the inhomogeneous solutions are computed by convolving them with the Teukolsky source. From the inhomogeneous solutions we can compute the radiated fluxes per mode from the complex asymptotic amplitudes, $C_{\ell m}^{\pm}$, of the radial Teukolsky solutions via Eq. (11).

Both in the source and in the orbital dynamics we linearize with respect to $\sigma$. Our resulting fluxes though are not linear in $\sigma$. This is because the radial Teukolsky equation contains a term proportional to the square of the mode frequency, $\omega=m\Omega$. In principle it would be possible to solve this equation to linear order in $\sigma$ but we have not attempted to do so in this work. Instead, at a range of fixed orbital frequencies, we numerically compute the fluxes for various values of $\sigma$ and fit this data using a polynomial and extract the linear in $\sigma$ piece. It is important to make this step as the quadratic and higher contributions to the raw frequency-domain flux data are not complete as we are not including, e.g., higher order corrections to the orbital dynamics.

The code is written in C++ and internally (particularly for the MST part of the calculation) it uses extended precision. We have high confidence in the results for three reasons: i) when $\sigma=0$ we recover known flux results non-spinning bodies to $\sim$14 significant digits ii) for $\sigma\neq 0$ and after extracting the linear in $\sigma$ contribution we see good agreement with the PN results – see Fig. 1 – and iii) we have reconstructed the metric perturbation at the particle using the standard CCK procedure Chrzanowski (1975); Cohen and Kegeles (1974); Kegeles and Cohen (1979) and from this we have shown that the local self-force experienced by the particle is balanced (to a relative error of $10^{-9}$) by the radiated fluxes through the spacetime boundaries. Details on each of these checks will be presented in Akcay et al. .

II.2.2 Numerical results: time-domain approach with Hamiltonian dynamics

The time-domain method solves the 2+1 Teukolsky equation obtained by separating Eq. (5) in the azimuthal ($\phi$) direction Harms et al. (2016b, a). The resulting wave equations for each $m$-mode are written in the scri-fixing hyperboloidal and horizon penetrating co-ordinates developed in Zenginoğlu (2011a, b); Harms et al. (2016a). The source of Teukolsky equation is implemented in a general way and then specified to Eq. (II.1.2) for a spinning test body (in the pole-dipole approximation) moving on an arbitrary trajectory Harms et al. (2016b); Zelenka et al. (2019). The Teukolsky equation is numerically solved using the method-of-lines with a 4th order accurate Runge-Kutta time integrator and 6th order accurate finite-differencing operators. The delta functions in the source can be discretized either using a narrow Gaussian or discrete delta functions; the former option is used in case of circular orbits as discussed in Harms et al. (2016b). The code was extensively tested and delivers multipolar waveforms and GW fluxes at null infinity with an accuracy well below the $1\%$ level up to $m=4$ modes  Nagar et al. (2014); Harms et al. (2016b, a); Lukes-Gerakopoulos et al. (2017). The data used for this work are produced exactly as described in Harms et al. (2016b, a). Circular equatorial trajectories for the test body are computed using the Hamiltonian formalism detailed in Sec. III of Harms et al. (2016a). The Hamiltonian of a spinning particle at linear order in the spin was originally obtained in Ref. Barausse et al. (2009). The same Hamiltonian can be recasted in certain specific effective-one-body (EOB) coordinates as illustrated in Ref. Bini et al. (2015). The EOB coordinates are transformed to the hyperboloidal ones via a transformation linear in $\sigma$. The EOB dynamics is compatible with the MPD dynamics with Tulczyjew and the Pirani Pirani (1956) SSCs across almost the whole spin and frequency range for the Schwarzschild background Harms et al. (2016a). One must, however, be aware that that the Hamiltonian circular dynamics and the MPD dynamics described above are fully equivalent and compatible only when the spin of the particle is small. We address the reader to Ref. Harms et al. (2016a) for additional details.

Following previous work of some of us Harms et al. (2016a); Lukes-Gerakopoulos et al. (2017), we want to accurately compute numerical fluxes for values of the spins that are large, i.e., $\sigma\sim 1$ and for orbits that are near the spin-dependent last-stable orbit. Past work Harms et al. (2016a) suggests that this can be done using the Hamiltonian formalism. Like the FD results mentioned above, our results will not be linear in $\sigma$, essentially for the same reasons mentioned above. Thus, to provide a consistent comparison with FD results, we also extracted the linear-in-sigma piece out of the numerical results.

II.3 Comparing numerical waveform amplitudes

Since the main aim of the paper is to check resummation procedure for gravitational waveform amplitudes, it is convenient to directly use these quantities for comparing the results obtained with the two numerical approaches. In this respect, the extraction of the linear piece mentioned in the previous section is now performed at the level of the waveform amplitudes (defined below) and not on the fluxes. Likewise for the fluxes, each waveform multipole is factorized in the product of a Newtonian (leading-order) contribution and a relativistic correction

[TABLE]

where ${\cal R}h_{\ell m}^{(N,\epsilon)}$ is the Newtonian prefactor that is analytically known (see Appendix B.1). In practice, it is convenient to focus only on the relativistic correction, since it is a function of order unity whose PN expansion has the structure $\hat{h}_{\ell m}^{(\epsilon)}\simeq 1+x+x^{3/2}+\dots$. Figure 2 offers a comprehensive comparison of the various numerical data at our disposal for an illustrative sample of multipoles and spins. The top part of each panel of the figure shows 4 curves: (i) the outcome of the time-domain code with Hamiltonian dynamics; (ii) the same quantity where one has subtracted the nonlinar-in-spin part; (iii) the outcome of the frequency-domain code with MPD dynamics: (iv) the same quantity where one has subtracted the nonlinear-in-spin part. The vertical line in the plot marks the location of the $\sigma$-dependent frequency of the last-stable-orbit (LSO), as obtained, for example, in Ref. Harms et al. (2016a)

[TABLE]

As can be seen from Fig. 2, the linear-in-spin results from the two codes are in the excellent agreement reaching a fractional accuracy of $\sim 10^{-5}$ (this is limited by the precision of the TD code and the need to extract the linear-in-spin contribution). Note that we do not expect agreement between the nonlinear in spin results as each code effectively includes different pieces of the nonlinear contribution (either through their dynamics or other aspects of the calculation). We emphasize that, to the best of our knowledge, this is the first successful comparison between two completely independent Teukolsky codes with a spinning secondary object. The excellent agreement between the two codes (after linearizing in spin) and the agreement between the FD and PN results presented in Fig. 1 gives us a high degree of confidence in our results.

Since the PN results are linear in the spin, we will use the linearized FD data as the target points to verify the accuracy of the various analytical approximations to the waveform amplitudes.

III Comparing analytical and numerical results

We turn now to comparing the “exact” numerical data computed in the previous section to different analytical approximations of $|\hat{h}_{\ell m}^{(\epsilon)}|$. More precisely:

• (i)

In Sec. III.1, we consider the straightforward PN-expanded expressions of $|\hat{h}_{\ell m}^{(\epsilon)}|$ up to 5.5PN order.

• (ii)

In Sec. III.2 we consider the, now standard, factorization and resummation of Refs. Damour et al. (2009); Pan et al. (2011a). Note that this analytical approach is adopted in state-of-the-art EOB waveform models for coalescing binaries Bohé et al. (2017); Nagar et al. (2018).

• (iii)

In Sec. III.3 we explore the performance of a hybrid approach that consists of: implementing the orbital factorization, and subsequent resummation, of Refs. Nagar and Shah (2016); Messina et al. (2018) for the $\epsilon=0$ modes, while the $\epsilon=1$ only benefit the factorization of the tail contribution.

• (iv)

In Sec. III.4 we test the special orbital-factorization and resummation of odd-$m$ modes that is adopted in the $\nu\neq 0$ case and that was recently shown to yield excellent consistency between EOB and state-of-the-art NR waveform amplitudes.

III.1 PN-expanded waveform amplitudes

The comparison with the plain PN-expanded waveform amplitudes is exhibited in Fig. 3. One observes good consistency between the numerical and analytical results at low frequencies, as expected. This progressively worsens towards the LSO, that now is indicated by colored markers on the figure. This worsening is prominent for higher modes. It is, however, interesting to note that the PN-expanded $m=1$ modes already offer an excellent representation of the numerical data up to the LSO and for any value of $\ell$.

III.2 Standard resummation: factorized and resummed amplitudes

We start by investigating whether the standard factorization and resummation procedure of waveform amplitudes of Refs. Damour et al. (2009); Pan et al. (2011b) can reduce the gap between the PN-expanded and numerical amplitudes towards the LSO seen in Fig. 3. The PN-expanded circularized waveform amplitudes are resummed in the following factorized form

[TABLE]

where $\epsilon$ is the parity of $\ell+m$, $\hat{S}^{(\epsilon)}$ is the $\mu$-normalized source of the field, $|\hat{h}_{\ell m}(x)|$ the modulus of the tail factor and $\rho_{\ell m}$ the residual factorized amplitudes. The squared modulus of the tail factor is given by Eq. (59) of Ref. Damour et al. (2009), and reads

[TABLE]

The source $S^{(\epsilon)}$ is either the $\mu$-normalized energy along circular orbits ($\epsilon=0$), or the Newton-normalized angular momentum ($\epsilon=1$). From the PN-expanded Newton-normalized fluxes $\hat{F}_{\ell m}$ we then calculate the PN-expanded residual amplitudes $\rho_{\ell m}$. To do so, we need the energy and angular momentum of a spinning body along circular orbits of Schwarzschild spacetime at linear order in the spin. These were obtained in Eqs. (81) and (82) of Ref. Harms et al. (2016a),

[TABLE]

where the expression for $\hat{j}$ is obtained multiplying by $\sqrt{x}$ (i.e., the inverse of the Newtonian angular momentum along circular orbits) Eq. (81) of Ref. Harms et al. (2016a). The PN-expanded $\rho_{\ell m}$’s are then given by

[TABLE]

where $T_{n}$ denotes a Taylor-expansion of order $n$ and $\hat{F}_{\ell m}$ are the Newton-normalized PN-expanded energy fluxes. The $\rho_{\ell m}$ functions are here written as the sum of an orbital (spin-independent) and a spin-dependent term

[TABLE]

The PN knowledge of $\rho_{\ell m}^{\sigma}$ before this work was limited to global 3.5PN (NNLO) and the functions are written explicitly in Ref. Messina et al. (2018); Cotesta et al. (2018). Here the computation of each $\rho_{\ell m}^{\sigma}$ is pushed up to relative 5.5PN order, i.e., next-to-next-to-next-to-next-to-leading-order, N4LO in the spin-orbit coupling. We list the functions explicitly in Appendix C. One verifies that the 3.5PN-accurate truncation of our results agrees in full with the corresponding formulas of Refs. Messina et al. (2018); Cotesta et al. (2018).

The analytical $|\hat{h}_{\ell m}|$’s from Eq. (15) are compared with the corresponding numerical ones in Fig. 4. The effect of the standard resummation can be summarized as follows. First, the analytically resummed $\ell=m$ modes deliver a rather good approximation to the numerical functions. One can see that this remains true up to $\ell=m=6$. However, we see that the procedure gives a rather unreliable result (even for small values of $x$) for the $\ell+m=\text{odd}$ modes. This illustrates that the standard resummation approach introduced in Ref. Pan et al. (2011b) can become highly inaccurate in some corners of the black-hole binary parameter space. As such, it should be replaced by something else that is more robust. One possibility is proposed in the next section.

III.3 Improved resummation

References Nagar and Shah (2016); Messina et al. (2018) presented an alternative factorization and resummation procedure based on the idea of first factoring out the orbital, spin-independent, contribution to the $\rho_{\ell m}$’s and then independently resumming, in various ways, the orbital and spin factors. The same procedure was implemented in similar ways for the odd-parity and even-parity modes. In the case of a nonspinning particle on circular orbits around a Kerr black hole, Refs. Messina et al. (2018) illustrated that this procedure yields a remarkable analytical/numerical agreement between the $\rho_{\ell m}$’s (and fluxes) up to the LSO also for extremal values of the spin of the black hole (see e.g. Figs. 1 and 2) of Messina et al. (2018)). In this Section we test this procedure on the even-parity modes. By contrast, we apply a different factorization, that only concerns the tail factor, to the odd-parity modes. Some of the resummation procedures of Refs. Messina et al. (2018) applied to some of the odd-parity modes are instead discussed in Sec. III.4 below.

III.3.1 $\ell+m$=even: orbital factorization

We implement the orbital factorization of the $\rho_{\ell m}$’s of Messina et al. (2018) to all even-parity multipoles up to $\ell=6$ included. Consistently with Messina et al. (2018), the orbital factors are taken at 6PN order111Though the $(3,1)$ mode is taken at 5PN only. and Padé resummed according to Table I of Messina et al. (2018). More precisely, we use the standard $P_{2}^{4}$ Padé for all the multipoles except the $(3,1)$ one, that is kept at 5PN with a $P_{2}^{3}$ Padé (the bold values in the table are not considerated). The only exception is the $(4,2)$ multipole, that in this paper is resummed with a $P_{2}^{4}$ Padé instead of a $P_{0}^{6}$ Taylor Series, since the difference between the two choices in this case is minimal. By resumming the spin-dependent factors by taking their inverse-Taylor representation, the resummed residual amplitudes finally read

[TABLE]

where $\hat{\rho}_{\ell m}^{\sigma}=T_{n}\left[1+\rho_{\ell m}^{\sigma}/\rho_{\ell m}^{\rm orb}\right]$ and we defined $\overline{\hat{\rho}_{\ell m}^{\sigma}}\equiv(T_{n}[(\hat{\rho}^{\sigma}_{\ell m})^{-1}])^{-1}$. Finally, the even-parity waveform amplitudes read

[TABLE]

The analytical/numerical agreement, displayed in Fig. 5, is essentially comparable to the standard approach shown in Fig. 4 above.

III.3.2 $\ell+m=$odd: factorizing the tail factor only

For the odd-parity modes we suggest here to follow a new route: factor out only the tail factor spin part of the modulus, $|h_{\ell m}^{\sigma}|$, while keeping the orbital part $|h_{\ell m}^{\rm orb}|$ factorized as usual (19), i.e., with the orbital angular momentum $\hat{j}^{\rm orb}$ factored out. Note that for this particular calculation, we keep $\rho_{\ell m}^{\rm orb}$ at 5PN accuracy and in Taylor-expanded form. The rationale behind the choice of not factorizing the orbital angular momentum is that, in the presence of a spinning body, the source of the field is given by the sum of two separate pieces, one proportional to $\hat{j}^{\rm orb}$ and another one to $\sigma$. This is in particular the structure of the source of the Regge-Wheeler-Zerilli equation for a spinning test-body. It seems then less sound to factor out $\hat{j}^{\rm orb}$ also from the $\sigma$-dependent term. Starting from

[TABLE]

we factorize each term separately as

[TABLE]

The resummed odd-parity waveform amplitudes finally read

[TABLE]

It is important to stress that, for simplicity, we are here not using Padé approximants of the orbital $\rho_{\ell m}^{\rm orb}$ in Eq. (26). The analytical and numerical amplitudes are compared in Fig. 5. The agreement between the two is remarkable, and way better than the one obtained with the standard approach shown in Fig. 4 above. To better quantify the agreement, we list in Table 1 the fractional differences computed at $x_{\rm LSO}$. Although not shown in the table, we have also verified that an analogous quantitative agreement holds for $\sigma=\pm 0.99$,

III.4 Resumming the dominant $\mathbf{m}$=odd modes

consistently with the comparable-mass case.

Up to this point we have seen that there are different procedures for resumming the waveform amplitudes depending on the parity of $\ell+m$. Actually, detailed studies of the waveform amplitudes for two objects of masses $(M_{1},M_{2})$ Nagar and Shah (2016); Messina et al. (2018) have illustrated that one should also carefully separate the analytic treatment depending on the parity of $m$ only. In particular, Refs. Nagar and Shah (2016); Messina et al. (2018) introduced a special analytical treatment of the orbital-factorized spin-dependent functions $\hat{f}_{\ell m}^{S}$ when $\nu\equiv M_{1}M_{2}/(M_{1}+M_{2})^{2}\neq 0$. More precisely, one shows that, when $\nu\neq 0$, these functions are naturally written as the sum of two separate Taylor expansions, one proportional to $X_{12}=X_{1}-X_{2}=\sqrt{1-4\nu}$ and the other to $\tilde{a}_{12}\equiv X1\chi_{1}-X_{2}\chi_{2}$, where $X_{i}\equiv M_{i}/(M_{1}+M_{2})$ and $\chi_{i}\equiv S_{i}/M_{i}^{2}$ are the dimensionless spins of the two objects. Each Taylor expansion can eventually be resummed taking its inverse-Taylor representation. This resummation of the $m$-odd waveform amplitudes has been recently incorporated in a new, multipolar, EOB waveform model. One of the remarkable features of this analytical choice is that the zero in the $(2,1)$ time-domain amplitude, that exists in certain regions of the parameter space, is quantitatively consistent with a similar feature found in state-of-the-art NR simulations. On top of this, Ref. Nagar and Shah (2016), pointed out that the resummation procedure is reliable and accurate also in the special case of a nonspinning particle orbiting a spinning black hole.

In this section we investigate to which extent the $m$-odd multipoles resummation approach remains robust and stable also in the case of a spinning test-body on a Schwarzschild black hole. To do so, we focus on the same $m$-odd modes considered in the $\nu\neq 0$ case, $(2,1)$, $(3,3)$, $(3,1)$, $(4,3)$, $(4,1)$. In addition, each mode is truncated at the largest PN order that carries known $\nu$-dependent corrections 222Even if we have more spinning-particle terms, we cannot use them since we cannot consistently split each PN order in the two separate series.. Consistently with Ref. Riemenschneider et al. , the $m$-odd waveform amplitudes are written as

[TABLE]

The Padé approximants of the orbital part $P^{m}_{n}[\rho_{\ell m}^{\rm orb}]$ are the same, adopted in Ref. Messina et al. (2018), that is $P^{5}_{1}$ for $(2,1)$, $P^{3}_{2}$ for $(3,1)$, $P^{6}_{0}$ for $(5,5)$ and $P^{4}_{2}$ for $(3,2)$, $(4,3)$ and $(4,1)$. On the other hand, the spin-dependent terms $\hat{f}_{\ell m}^{\sigma}$ are given by

[TABLE]

Following Ref. Messina et al. (2018), the overbar indicates that each term is resummed using its inverse-Taylor representation. The function written above explicitly read

[TABLE]

The analytical $|\hat{h}_{\ell m}|$ are compared to the numerical ones in Fig. 6. The most interesting result displayed in the figure is that both the $(2,1)$ and $(3,3)$ analytical amplitudes deliver a reasonably accurate representation of the numerical data up to the LSO, although this is not as good as the tail-factorized case discussed in Fig. 5. Note however that in that case we were using spin information truncated at 5.5PN accuracy. Thus, to produce a meaningful comparison we need to redo the calculation of Fig. 5 truncating the $\ell=2$, $m=1$ mode at 2.5PN. By contrast, the orbital $\rho_{21}^{\rm orb}$ function is kept at 5PN accuracy and in Taylor-expanded form. This new comparison is displayed in Fig. 7. Interestingly, despite the reduce PN information, the analytical/numerical agreement is visibly better than the one displayed in the top-left panel of Fig. 6. We also performed the same analysis for the $\ell=4$ modes. Although we found an improvement, the truncation at 1.5PN of the spin information is not sufficient to provide a good agreement up to the LSO location.

In conclusion, the result of Fig. 7 gives further support to the need of exploring the performance of the tail-only factorization also for comparable mass binaries.

IV Conclusions

In this paper we have collected several new results concerning the gravitational waveform fluxes and amplitudes emitted by a spinning particle (i.e., a spinning test-black hole) on circular orbits around a Schwarzschild black hole at linear order in the particle spin. Our main findings are summarized as follows:

• (i)

We have analytically computed the PN-expanded multipolar energy fluxes up to $\ell=7$. Each multipole is obtained at 5.5 PN order beyond the leading, Newtonian, contribution. This improves our current analytical knowledge of the fluxes in this corner of the parameter space, that was previously known only up to global 3.5PN order.

• (ii)

We have computed numerically the multipolar energy fluxes (and waveform amplitudes) up to $\ell=6$ using two different, and independent, numerical codes. One is a frequency-domain code; the other uses a time-domain approach. We demonstrated the excellent mutual consistency between the two numerical approaches once the numerical data are suitably linearized in the spin of the particle. This allows us to provide accurate, circularized, waveform amplitudes up to the LSO for different values of the particle spin.

• (iii)

We compared to these numerical data, considered as exact, the PN-expanded analytical amplitudes as well as different flavors of their resummation. In particular, we could show that the standard resummation approach of Ref. Pan et al. (2011b), that is implemented in state-of-the-art waveform models for coalescing black-hole binaries, such as SEOBNRv4 Bohé et al. (2017) or TEOBResumS Nagar et al. (2018), is inaccurate for $\ell+m$ odd modes and provides nonnegligible differences with the numerical data already in the early inspiral. This indicates that the procedure of Ref. Pan et al. (2011b) should be improved to construct a waveform model robust and accurate all over the BBH parameter space.

• (iv)

The most important result of this work is that the factorization of only the tail factor in the $\ell+m=\text{odd}$ modes allows one to obtain an excellent analytical/numerical agreement up to the LSO (and even below), for $-0.9\leq\sigma\leq+0.9$, and for all odd-parity modes analyzed up to $\ell=5$, see Fig. 5. This result opens at least two questions. First of all, one wonders whether the only-tail factorization would be helpful also to resumming the fluxes emitted by a nonspinning particle orbiting a Kerr black hole. Resumming that PN series has always been a challenge Pan et al. (2011b) that was solved, to a certain extent, only through the orbital factorization (and additional resummation) procedure of Refs. Nagar and Shah (2016); Messina et al. (2018). Hopefully the only-tail factorization approach might be useful in that case at least for what concerns the terms linear in the black-hole spin. On the other hand, our findings also suggest that it might be worth investigating the performance of the only-tail factorization procedure also in the $\nu\neq 0$ case. The accuracy of such odd-parity waveform amplitudes should then be carefully evaluated by performing extensive comparisons between a so-constructed EOB waveform model and state-of-the-art NR simulations.

• (v)

We also explored the accuracy of the orbital-factorization/resummation approach of Ref. Nagar and Shah (2016); Messina et al. (2018), notably in the form developed for the odd-$m$ modes in the $\nu\neq 0$ case. We do so by truncating the PN information at the same PN order where $\nu$-dependent corrections are known. This approach has been recently used to improve the behavior of the $\ell=2$, $m=1$ waveform amplitude in a new, EOB-based, multipolar waveform model Riemenschneider et al. . For this mode, we found (see top-left panel of Fig. 6) that the numerical/analytical agreement is rather good, though not at the level of the only-tail factorized case mentioned above. The situation is even worse for the subdominant modes (except the $\ell=m=3$ one). Our results seem to suggest that the resummation of Ref. Nagar and Shah (2016); Messina et al. (2018) should be replaced by the only-tail factorization one, since it is simpler and more accurate. More investigations, notably when both objects are spinning, will be needed to confirm this preliminary conclusion.

Acknowledgements.

FM thanks IHES for hospitality at various stages during the development of this work. GLG is supported by Grant No. GAČR-17-06962Y of the Czech Science Foundation that is gratefully acknowledged. NW gratefully acknowledges support from a Royal Society - Science Foundation Ireland University Research Fellowship. SB acknowledges support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. We acknowledge the networking support received by the COST Action CA16104.

Appendix A Interpreting the spinning-body limit using the EOB dynamics

In the main text we have considered $-1\leq\sigma\leq 1$. Due to the definition of $\sigma$ in Eq. (3) we can also write

[TABLE]

This parameter, per se, can be thus meaningful also for the equal-mass case, with the understanding that the MPD approach (or the Hamiltonian approach) plus the solution of the Teukolsky equations for the fluxes can only deliver some of the contributions at leading-order in the mass ratio. More precisely, the value $\sigma=\pm 1$ can be interpreted as the case of an extremally spinning black hole, $s/\mu^{2}=\pm 1$ in an equal-mass binary, $\mu/M=1$. Although in this case the perturbative approach we are using is expected to be meaningless, in practice it is useful to compute, and test numerically, some of the terms that enter the complete PN expansion.

The euristic comparable-mass limit suggested by Eq. (43) can be put on a more solid ground starting from the complete EOB Hamiltonian for a two-body system with masses $(m_{1},m_{2})$, with the convention that $m_{1}\geq m_{2}$. We also refer the reader to Sec. III of Ref. Harms et al. (2016a) for complementary information. The complete EOB Hamiltonian reads

[TABLE]

where now it is $\mu\equiv m_{1}m_{2}/(m_{1}+m_{2})$, and $\nu=\mu/M$ is the symmetric mass ratio and

[TABLE]

is the effective Hamiltonian. The spin-orbit sector is given by

[TABLE]

where $p_{\varphi}=P_{\varphi}/\mu$ is the dimensionless orbital angular momentum, we introduced the symmetric spin combinations

[TABLE]

and $(G_{S},G_{S_{*}})$ are the gyro-gravitomagnetic functions. For the case where only the secondary, $m_{2}$ is spinning ($S_{1}=0$), $\hat{H}_{\rm SO}^{\rm eff}$ becomes

[TABLE]

After introducing powers of $m_{2}$ so to explicitly have a dimensionless spin variable $\chi_{2}\equiv S_{2}/m_{2}^{2}$ one obtains

[TABLE]

The extreme-mass-ratio limit is now defined by the condition $\nu\to 0$. In this limit, one sees that the term proportionalt to $G_{S}$ is suppressed with respect to the one proportional to $G_{S_{*}}$ because of the different $\nu$ dependence. In the $\nu\to 0$ limit one finally has that the spin-orbit sector of the real spin-orbit Hamiltonian describing a spinning test-body orbiting a Schwarzschild black hole is given by

[TABLE]

where $G_{S_{*}}$ is known analytically in closed, non-PN expanded, form Barausse et al. (2009); Barausse and Buonanno (2010); Bini et al. (2015). One then defines the dimensionless spin $\sigma$ as

[TABLE]

Thus, since $-1\leq\chi_{2}\leq+1$, and $0\leq\nu\leq 1/4$, one concludes, from the EOB (and thus PN) perspective, that $-1/4\leq\sigma\leq 1/4$. The usual limit $-1\leq\sigma\leq 1$ is obtained writing $\sigma$ as

[TABLE]

and then addtionally expanding the symmetric mass ratio in powers of the mass ratio $m_{2}/m_{1}\approx\mu/M\ll 1$.

Appendix B Spinning particle on Schwarzschild background: PN-expanded multipolar energy fluxes

B.1 Newtonian prefactors

Each multipolar flux is written as the product of a Newtonian (or leading-order) prefactor and of a PN correction, i.e.

[TABLE]

where $\epsilon=0,1$ indicates the parity of $\ell+m$. For each multipole, the Newtonian prefactor can be written in closed form as

[TABLE]

where the Newtonian waveform ${\cal R}h_{\ell m}^{(N,\epsilon)}$ is explicitly given as Damour et al. (2009)

[TABLE]

with

[TABLE]

and

[TABLE]

B.2 Newton-normalized PN fluxes

The PN-expanded fluxes up to $\ell=7$ are presented below. For completeness and future reference we also keep all the nonspinning terms. All of the series given below can be found digitally in the PostNewtonianSelfForce package in the Black Hole Perturbation Toolkit BHP .

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Multipolar 5.5PN EOB relativistic residual amplitudes derived from the new spinning particle on Schwartzschild fluxes results

The spin-dependent part of the PN-expanded residual relativistic amplitudes obtained from the fluxes of Sec. B.2 read

[TABLE]

[TABLE]

[TABLE]

[TABLE]