Eccentric Gravitational Wave Bursts in the Post-Newtonian Formalism

TL;DR

This paper develops a detailed post-Newtonian framework to model eccentric gravitational wave bursts, enhancing detection strategies for highly eccentric binaries in gravitational wave astronomy.

Contribution

It introduces a generic post-Newtonian formalism for constructing time-frequency priors for eccentric GW bursts, including a third post-Newtonian order model, the most accurate analytic model to date.

Findings

Provides a new formalism for eccentric GW burst modeling

Generates a third post-Newtonian order burst model

Enhances detection strategies for eccentric binaries

Abstract

The detection of GW150914 by ground based gravitational wave observatories has brought about a new era in astrophysics. At optimal sensitivity, these observatories are expected to detect several events each year, with one or two of these occurring with non-negligible eccentricity. Such eccentric binaries will emit bursts of gravitational radiation during every pericenter passage, where orbital velocities can reach greater than ten percent the speed of light. As a result, such binaries may prove to be powerful probes of extreme gravitational physics and astrophysics. A promising method of achieving detection of such binaries is through power stacking, where the power in each burst is added up in time-frequency space. This detection strategy requires a theoretical prior of where the bursts will occur in time and frequency so that one knows where to search for successive bursts. We here present a generic post-Newtonian formalism for constructing such time-frequency model priors at generic post-Newtonian order. We apply our formalism to generate a burst model at third post-Newtonian order, making it potentially the most accurate, fully analytic model to date.

Tables

Table 1. Table 1: Initial values of the PN expansion parameter x 𝑥 x for the set of compact binary systems studied. The values are obtained by requiring the initial GW frequency to be 10 ​ H ​ z 10 H z 10{\rm Hz} . The final column provides an estimate of the semi-major axis of the binary, since a r = M / x + 𝒪 ​ ( 1 ) subscript 𝑎 𝑟 𝑀 𝑥 𝒪 1 a_{r}=M/x+{\cal{O}}(1) in PN theory.

System	$m_1\left[M_{\odot }\right]$	$m_2\left[M_{\odot }\right]$	$e_{t,0}$	$x_0$	$1/x_0$
NSNS	$1.4$	$1.4$	$0.9$	$7.35\times {10}^{-4}$	$1360$
NSBH	$1.4$	$10$	$0.9$	$1.85\times {10}^{-3}$	$541$
BHBH	$10$	$10$	$0.9$	$2.67\times {10}^{-3}$	$375$

Equations

$$\frac{(t_i-t_{i-1}{)}_{\mathrm{PN}}}{(t_i-t_{i-1}{)}_{\mathrm{N}}}$$

$$\frac{f_i^{\mathrm{PN}}}{f_i^{\mathrm{N}}}$$

$$\frac{\delta t_i^{\mathrm{PN}}}{\delta t_i^{k\mathrm{N}}}$$

$$\frac{\delta f_i^{\mathrm{PN}}}{\delta f_i^{\mathrm{N}}}$$

$${\left[1+\vec{R}(e_t,\eta ;v_p)\cdot \vec{X}(v_p)\right]}^{(-1)}\doteq 1+{\vec{R}}^{(-1)}(e_t,\eta ;v_p)\cdot \vec{X}(v_p),$$

$$\frac{(v_{p,i}-v_{p,i-1}{)}_{\mathrm{PN}}}{(v_{p,i}-v_{p,i-1}{)}_{\mathrm{N}}}$$

$$\frac{(\delta e_{t,i}-\delta e_{t,i-1}{)}_{\mathrm{PN}}}{(\delta e_{t,i}-\delta e_{t,i-1}{)}_{\mathrm{N}}}$$

$$E$$

$$L$$

$$E_i$$

$$L_i$$

$$\Delta E_{(i,i-1)}$$

$$\Delta L_{(i,i-1)}$$

$$\Delta E_{(i,i-1)}$$

$$\Delta L_{(i,i-1)}$$

$$E_i$$

$$L_i$$

$$r_{p,i}$$

$$\delta e_i$$

$$P=\frac{2\pi r_p^{3/2}}{M^{1/2}(1-e{)}^{3/2}}.$$

$$t_i=t_{i-1}+\frac{2\pi }{M^{1/2}}{\left[\frac{r_{p,i}\left(r_{p,i-1},\delta e_{i-1}\right)}{\delta e_i\left(r_{p,i-1},\delta e_{i-1}\right)}\right]}^{3/2}$$

$$t_i-t_{i-1}={\int }_0^{2\pi }\frac{d\psi }{\dot{\psi }}.$$

$$f_{\mathrm{GW}}=\frac{1}{2\pi {\tau }_{\mathrm{GW}}},$$

$${\tau }_{\mathrm{GW}}$$

$$v_p=\sqrt{\frac{M(1+e)}{r_p}}$$

$${\tau }_{\mathrm{GW}}$$

$$f_i=\frac{M^{1/2}{\left[2-\delta e_i\left(r_{p,i-1},\delta e_{i-1}\right)\right]}^{1/2}}{2\pi {\left[r_{p,i}\left(r_{p,i-1},\delta e_{i-1}\right)\right]}^{3/2}}.$$

$$\delta t={\xi }_t{\tau }_{\mathrm{GW}}$$

$$\delta f={\xi }_f{f}_{\mathrm{GW}}$$

$$\delta t_i$$

$$\delta f_i$$

$$v_{p,i}$$

$$\delta e_{t,i}$$

$$\Delta v_{p,(i,i-1)}$$

$$\Delta e_{t,(i,i-1)}$$

$$v_{p,i}$$

$$\delta e_{t,i}$$

$$t_i$$

Full text

1 Introduction

Over the past several years, studies have shown that exotic formation channels could lead to a population of highly eccentric compact binaries whose gravitational wave (GW) emission would be in the band of current ground-based detectors [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. One such formation channel is dynamical captures in dense stellar environments, such as globular clusters and galactic nuclei [1, 2, 10]. In such dense environments, compact objects initially on hyperbolic trajectories can become bound after passing through closest approach due to GW emission or tidal interactions, with the subsequent bound binary having high eccentricity ($e\sim 1$). On the other hand, the Kozai-Lidov mechanism [3, 11, 4, 5, 6, 7, 8, 9, 12], and other three-body interactions [7] in hierarchical triple systems, can induce resonances that drive the inner binary to the parabolic limit.

Although expected to be rare [13], eccentric binaries could prove to be powerful probes of astrophysical dynamics. Event rates of eccentric inspirals due to dynamical captures have wide error bars, typically around two orders of magnitude [1]. Such large error is largely due to the unconstrained populations of compact objects within dense environments [13]. Likewise, the tightening of binaries due to three-body interactions in galactic nuclei are similarly uncertain [14]. Hence, the detection of GWs from eccentric inspirals would provide information about the mass function of black holes (BHs) and neutron stars (NSs) in these environments, allowing us to probe astrophysics that has proven difficult to extract from electromagnetic observations [15, 16, 17, 18].

Another promising area of interest for eccentric binaries is testing Einstein’s theory of General Relativity (GR). For highly eccentric binaries, the distance of closest approach can be small relative to the semi-major axis of the orbit, leading to systems with pericenter velocities greater than $10\%$ the speed of light. At such high velocities, the GW luminosity in each burst will typically be $\sim(10^{-4}-10^{-3})L_{{\mbox{\tiny Pl}}}$, where $L_{{\mbox{\tiny Pl}}}$ is the Planck luminosity. For comparison, a BH-BH, quasi-circular binary emits radiation at $\lesssim 10^{-5}L_{{\mbox{\tiny Pl}}}$ during the early phase of inspiral, increasing rapidly close to merger and eventually reaching $\sim 10^{-2}L_{{\mbox{\tiny Pl}}}$ only at merger. As such, the GWs from eccentric binaries could capture effects from the extreme gravity regime (i.e. where gravity is dynamical, strong and non-linear) during many pericenter passages [19].

These GWs may be detected by current and upcoming ground-based GW detectors, and hence the study of GWs from eccentric binaries has never been more urgent. The advanced Laser Interferometer Gravitational Wave Observatory (aLIGO) [20] has already achieved the first detection of GWs with the event GW150914 [21]. The advanced VIRGO Interferometer (aVIRGO) [22] will be coming online in 2016-2017, with additional detections expected during this period. Plausible estimates of the event rates of the inspiral of compact object binaries predict that these detectors could see $\sim 5-10$ events per year in the very near future [13, 23, 24]. Based on our current knowledge of the formation channels of eccentric compact binaries [1, 2, 3, 4, 5, 6, 7, 8, 7, 9], we might expect that one or two of these events will enter the LIGO band with non-negligible orbital eccentricity. Looking toward the future, KAGRA [25, 26] in Japan is expected to come online by the end of this decade and LIGO India [27] in the beginning of the next decade. Once these detectors are operational, the number of detected events per year will necessarily increase, thus increasing the probability of detecting eccentric inspirals.

The typical strategy for detecting well-understood GWs with ground-based interferometers is to use matched filtering [28, 29, 30, 31, 32, 33, 34]. Effectively, a set of templates that best describe the signal buried in the data are used to extract the latter and estimate its physical parameters. This detection strategy hinges on having very accurate templates, as a small dephasing between the signal and a template can result in complete loss of detection [35, 36, 37, 38]. Hence, we must have a prior model of what to search for in the detector data output. For highly eccentric binaries, the GW emission resembles a set of discrete bursts somewhat localized in time and frequency. These bursts are centered around pericenter passage, where the orbital velocity is highest and where the binary spends the least amount of time, and thus very little GW power is contained within each individual burst. This issue alone makes matched filtering a rather impractical search strategy for such GWs, but it is further compounded by the fact that there are few, fully-analytic and accurate templates for highly eccentric binaries with which to perform matched filtering computationally efficiently [39, 40].

An alternative search strategy was presented in [41] using a power stacking method. Ideally, if one could register a set of bursts in a time-frequency decomposition of the data stream, one could then stack the power within each burst, thus creating an enhanced data product. For $N$ bursts with the same signal-to-noise ratio (SNR), the amplification in the SNR relative to the SNR in a single burst would scale as $N^{1/4}$. Although sub-optimal compared to matched filtering, the power stacking method is more robust to modeling errors and more efficient in detection of eccentric signals than current un-modeled burst methods [41]. Power stacking, however, still requires a model of where the bursts will occur in time-frequency space given some initial starting point, with which to construct a prior to search for successive bursts.

Such a burst model was developed in [19] for tracking the bursts in time-frequency space. In general, a burst model is one that treats the bursts as N-dimensional objects in the detector’s data stream and tracks the geometric centroid and volume of the bursts from one to the next as the binary inspirals. To do this, [19] considered Keplerian orbits perturbed by quadrupolar gravitational radiation. We will refer to this model as Newtonian in the sense that it is obtained from a fully relativistic model expanded to lowest, non-vanishing order about small orbital velocities and weak gravitational fields. The benefit of working in such a simplified scenario was that it allowed for the exploration of whether such burst signals could be used to test well-motivated deviations from GR with eccentric signals as a proof-of-concept.

The ultimate goal, of course, is to create a model that is as accurate as possible relative to the signals present in nature. The modeling of the coalescence of compact objects in full GR is an exceedingly difficult problem, which has only been solved numerically (predominantly for quasicircular binaries) in the passed several years. For eccentric binaries, numerical simulations in full GR are more computationally expensive, and thus, to obtain only a few orbits at the desired numerical accuracy requires much more computational time than that needed in the evolution of quasicircular binaries. At least for the moment, the pure numerical modeling of eccentric binaries over the last thousand orbits in full GR is currently an intractable problem. On the other hand, we could work to extend the Newtonian burst model by considering relativistic corrections to Newtonian dynamics.

The post-Newtonian (PN) formalism [42, 43, 44, 45, 46, 47] allows for a systematic treatment of $v/c$ corrections to Newtonian dynamics, where $v$ is the orbital velocity and $c$ is the speed of light. For bound binaries, the orbital velocity is connected to the gravitational field strength through a Virial relation: a term of ${\cal{O}}(v^{2}/c^{2})$ is comparable to a term of ${\cal{O}}[GM/(Rc^{2})]$, where $M$ and $R$ are the characteristic mass and orbital separation of the system, and $G$ is Newton’s gravitational constant. Hence, the PN formalism for binaries is simultaneously a post-Minkowskian expansion, i.e. it is both an expansion in $v/c\ll 1$ and an expansion in $GM/(Rc^{2})\ll 1$.

The PN formalism has had a wide range of success in the modeling of binaries and their GW emission. At present, the GW emission from quasicircular binaries has been calculated to 3.5 PN order111A term proportional to $(v/c)^{2N}$ relative to its controlling factor will be said to be of N PN order. [46] and to 4PN order in the effective one-body Hamiltonian [48]. To leading order in the mass ratio, the radiation fluxes at spatial infinity are currently known to 22 PN order [49], due to the formulaic nature of the calculation. The PN corrections to Newtonian dynamics for eccentric binaries have proven more difficult to calculate. Damour and Deruelle [50, 51] found a Keplerian parameterization of the solution to the 1PN order equations of motion in terms of the eccentric anomaly $u$. This quasi-Keplerian (QK) representation was extended to 2PN order in [52] and to 3PN order in [53].

The QK parametrization must be enhanced to include dissipation if one wishes to obtain accurate waveform models. This parameterization is a purely conservative representation of the orbital motion of eccentric binaries, because the orbital energy and angular momentum are assumed to be conserved. Dissipation occurs because GWs carry energy and angular momentum away from the system, causing the binary to inspiral and eventually merge. For eccentric binaries, the GW energy and angular momentum fluxes have been computed to full 3PN order [54, 55]. With these fluxes at hand and assuming small eccentricities ($e\lesssim 0.1$), Refs. [56, 57, 58, 59] constructed time-domain and frequency domain waveforms to 2PN order. Waveform templates also exist for higher eccentricity systems ($e\lesssim 0.4$), for example through the hybrid time-domain x-model of [39] and the hybrid frequency-domain model of [40]. These models, however are really not applicable to highly eccentric binaries.

The burst model previously developed [19] is currently the only purely analytic model for highly elliptic orbits. We say the orbits are highly elliptic, and not highly eccentric because the latter implies that the eccentricity could be large and potentially greater than unity. On the other hand, “highly elliptic” indicates that we are always considering bound orbits. The burst model has currently only been developed to Newtonian order, which is an artifact of our desire to simplify our previous analysis as much as possible to be able to consider tests of GR.

Nature, however, is not Newtonian, and thus, extending the burst model into the PN formalism serves two purposes: to improve its potential of aiding in the detection of highly elliptic binaries and to enhance its ability to perform interesting and important science. The burst model was conceived with the idea of testing GR. However, if there is an inherent modeling error within the GR model, it is possible that such modeling errors could fool us into believing we have detected a non-GR signal if we are not careful [35, 36]. Furthermore, the power stacking method is not immune to modeling error in the detection of signals and the extraction of their parameters [41]. Hence, the detection of such binaries and extraction of important astrophysics hinges on having an accurate prior to predict where the bursts occur in time-frequency space.

1.1 Executive Summary

We here extend the burst model developed in [19] to higher PN order. The Newtonian burst model in [19] focused on the bursts emitted during the inspiral of the binary only. Similarly, we here also focus on the inspiral of highly elliptic binaries within the PN framework. The motivation for developing a purely analytic model of the inspiral is the potential for the later construction of phenomenological inspiral-merger-ringdown models, a quasi-circular version of which played a pivotal role in the first gravitational wave observations by aLIGO [21]. Here, we treat the bursts as two dimensional regions of excess power in a time-frequency decomposition of a detectors data stream. We treat the bursts as boxes with characteristic time and frequency widths, which allows for a discretization of the time-frequency decomposition into tiles, with the burst being those tiles that contain excess power [21, 41]. As with the Newtonian burst model, we characterize the sequence of bursts using the time and frequency centroids of the bursts, as well as the widths of the burst tiles, or alternatively the volume of the tiles, used to capture a certain amount of power within each bursts. These time-frequency observables are supplemented by a model describing the orbital evolution of the binary as a set of discrete, osculating Keplerian ellipses.

Similar to how the parameterized post-Einsteinian (ppE) burst sequence of [19] was a parameterized deformation of a simplified GR sequence, the PN burst sequence will be a parameterized deformation of the Newtonian order sequence. The deformations will scale with an increasing power of a particular PN expansion parameter, which we choose to be the pericenter velocity $v_{p}$. The coefficients of a $k/2$-PN order term, which scales as $v_{p}^{k}$, are then a set of functions $[{\cal{V}}_{k},{\cal{D}}_{k},P_{k},R_{k}]$ which are dependent on the physical parameters of the compact binary system. These functions correspond to the PN corrections to the rates of change of pericenter velocity and time eccentricity, and the expressions for the orbital period and pericenter distance, respectively. In this work, we neglect the spin of the compact objects and work in a point particle limit, such that these functions are only dependent on the time eccentricity $e_{t}$ and the symmetric mass ratio $\eta$. Hence, when working to $k/2$-PN order, one needs $4k$ functions $[{\cal{V}}_{k}(e_{t},\eta),{\cal{D}}_{k}(e_{t},\eta),P_{k}(e_{t},\eta),R_{k}(e_{t},\eta)]$ to parametrize all of the PN defomations.

We parametrize the PN burst sequence in time-frequency space by

[TABLE]

where $(t_{i},f_{i})$ are the centroid of the bursts and $(\delta t_{i},\delta f_{i})$ are the width and height of the tiles. We create the amplitude vector fields $[\vec{P},\vec{R}]$, whose components are the functions $[P_{k}(e_{t},\eta;v_{p}),R_{k}(e_{t},\eta;v_{p})]$, which we further specify as implicit functions of the PN expansion parameter $v_{p}$ since their form changes depending on the choice of expansion parameter. The components of $\vec{R}^{(-1)}$ are defined such that

[TABLE]

where the equality $\doteq$ should be understood as working in the limit of $v_{p}\ll 1$. The state vector $\vec{X}(v_{p})$ contains the powers of $v_{p}$ that characterized each PN order corrections, specifically $\vec{X}(v_{p})=(v_{p},v_{p}^{2},...,v_{p}^{k})$. Hence the dot products provide the complete sum of all terms in a PN expansion up to $k/2$-PN order.

These time-frequency burst parameters, specifically $(t_{i},f_{i},\delta t_{i},\delta f_{i})$, are functions of the symmetric mass ratio and the total mass of the binary, which are constant in time, as well as the pericenter velocity $v_{p,i}$ and eccentricity $e_{t,i}$ during each burst, which are evolving in time under the influence of radiation reaction. Hence, we must supplement the time-frequency sequence described above with the orbital evolution of the binary. To do this, we apply an osculating approximation that assumes the bursts are emitted instantaneously at pericenter, forcing the binary to move along a discrete set of Keplerian ellipses that osculate onto one another. The parameters of the $i$-th orbit will be functions of the parameters of the previous orbit, specifically

[TABLE]

where we have introduced the two new amplitude vector fields $[\vec{\cal{V}},\vec{\cal{D}}]$. In this work, we provide explicit expressions for the amplitude vector fields $\lambda_{\rm PNburst}^{a}\equiv(\vec{P},\vec{R},\vec{\cal{V}},\vec{\cal{D}})$ complete to 3PN order. Equations (1)-(4) and (6)-(7) provide the complete PN burst model, which we use to calculate the burst model to 3PN order using the results for $\lambda_{\rm PNburst}^{a}$.

How does this new burst model aid us in the detection of highly elliptic binaries? In a realistic search, the burst model acts as a prior on where the bursts will occur in time-frequency space. For example, once a search detects a burst of power within an interferometer data stream (even if this burst of power is not “loud” enough to allow to claim detection), the burst model can then be used to search over ”future” time-frequency space for successive bursts, as well as ”past” time-frequency space for bursts that may have been missed by previous searches. Physically, this amounts to searching over the parameters of the system that determine the prior, specifically the eccentricity and pericenter velocity during the initially detected burst, and the chirp mass and symmetric mass ratio of the binary.

The structure of the PN burst model should not be surprising given the structure of the ppE burst model in [19]. The ppE model requires four exponent parameters $a_{i}$, which govern the power of $v_{p}$ of the corrections, and four amplitude parameters $\alpha_{i}$, which depend on the coupling constants of the theory and the eccentricity of the binary. In the PN formalism, the exponent parameter $k$ becomes a known quantity and only changes when one goes to higher order in the expansion variable. The amplitude parameters have now been replaced with four amplitude vector fields222These are not true vector fields, but are a set of scalar functions that have been combined into a discrete sequence. The terminology used for these functions goes along with the notation we have used to simplify some of the expressions in this work. that parametrize the eccentricity and mass dependence of specific PN terms. As a result of this, rather than needing eight parameters as was the case in the ppE model, we require 4$k$ functions when working to $k/2$-PN order. These amplitude functions only depend on the initial eccentricity and the symmetric mass ratio, which together with the initial pericenter velocity and the chirp mass of the binary are the only parameters needed to define the model. Further, the fact that we require four vector fields to describe the burst model is a result of the fact that the model in only parametrized by four quantities: the orbital energy and angular momentum, and the energy and angular momentum fluxes of the GWs emitted by the system. Alternatively, as we will show, a different set of four parameters can be used: the orbital period, the mapping between the pericenter distance and pericenter velocity, and the rates of change of eccentricity and pericenter velocity due to radiation reaction. Working with these four quantities significantly improves the ease with which burst models can be constructed.

The remainder of this paper is dedicated to deriving the results presented above. In Section 2, we review the Newtonian order burst model and present a simplified formalism used for the construction of the PN burst model. Section 3 is dedicated to constructing a burst model at arbitrary PN order, which we later specialize to the cases of 1PN, 2PN, and 3PN orders. Section 5 discusses the results of the paper and their importance for future research. In this paper, we use geometric units where $G=c=1$.

2 Constructing Burst Models

This section is dedicated to reviewing how to create a burst model and the elements that go into such a model. We begin by reviewing the Newtonian burst model and how it was constructed. We then describe a new method of constructing burst models in general without any assumptions of the regime or theory of gravity we are working in. This new method greatly simplifies the construction of burst models, and will allow us to develop a completely generic GR PN model in the next section.

2.1 The Newtonian Burst Model

How do we construct a burst model? Recall that in Section 1 we defined a burst model as a theoretical model prior to describe how the bursts evolve in time-frequency space. Such a model would tell us how the centroid and the volume of the bursts evolve in time and frequency from one burst to the next. But this evolution depends on the orbital parameters of the system, which themselves are also changing in time due to dissipative effects, such as the emission of GWs. Therefore, a complete burst model must provide a one-to-one mapping between the evolution of the system’s physical parameters and how the bursts evolve in time and frequency. This requires the following ingredients:

1. I.

Orbital Evolution: A mapping that prescribes the evolution of the orbital parameters from one orbit to the next, including GW radiation-reaction. 2. II.

Centroid Mapping: A mapping that provides the time-frequency centroid of the burst $(t_{i},f_{i})$, given the centroid of the previous burst $(t_{i-1},f_{i-1})$, in terms of the orbital parameters of the system. 3. III.

Volume Mapping: A mapping that describes how the time-frequency volume of the bursts changes from one to the next, in terms of the orbital parameters of the system.

2.1.1 Orbit Evolution

Let us start by reviewing how the ingredients listed above can be computed to leading (i.e. Newtonian) order, focusing first on ingredient I (the orbital evolution). In Newtonian gravity, the orbital motion of two test particles can be described though Keplerian ellipses, which are characterized by two conserved quantities, the orbital energy $E$ and the orbital angular momentum $L$. Alternatively, one can parameterize any such orbit in terms of its pericenter distance $r_{p}$ and its orbital eccentricity $e$, which are related to $E$ and $L$ at Newtonian order in a PN expansion by

[TABLE]

where $\eta=m_{1}m_{2}/M^{2}$ is the symmetric mass ratio and $M=m_{1}+m_{2}$ is the total mass of the system.

The burst model requires knowledge of how $(E,L)$ or $(r_{p},e)$ evolve from one orbit to the next. Due to the nature of the emission of gravitational radiation in highly elliptic systems, we may treat the problem as a set of Keplerian orbits that change effectively instantaneously at pericenter, allowing the orbits to osculate onto one another. Hence, we may write

[TABLE]

where $\Delta E_{(i,i-1)}$ and $\Delta L_{(i,i-1)}$ are the changes in orbital energy and angular momentum due to GW emission from one orbit to the next, and the labels represent which orbit the above quantities are evaluated on. By “osculating orbits,” we mean that the elements of the Keplerian orbit are constant throughout the orbit except at pericenter, where they change drastically and the new elements define a new Keplerian orbit. In this approximation, one thus treats the radiation, and all changes generated by it, as arising instantaneously at pericenter.

In general, the total change in energy and angular momentum between times $T_{i-1}$ and $T_{i}$ due to GW emission is given by

[TABLE]

where $\psi$ is the true anomaly, $T_{i-1}$ and $T_{i}$ are the times of consecutive pericenter passages, and the dot refers to derivatives with respect to time. At Newtonian order, the GW energy and angular momentum fluxes, $\dot{E}$ and $\dot{L}$, are given, for example, by Eq. (12.78) in [47]. The fluxes are functions of the orbital elements, which are themselves functions of time through the true anomaly $\psi$. Thus, the above definitions would have to be supplemented with the time evolution of $\psi$ itself, namely $\dot{\psi}\left(r_{p},e,\psi\right)$. In addition, since the fluxes depend on the true anomaly, they contain gauge-dependent terms [47]. However, these terms vanish upon integration, leaving $\Delta E_{(i,i-1)}$ and $\Delta L_{(i,i-1)}$ independent of the radiation reaction gauge.

To evaluate the integrals in Eqs. (12) and (13), we perform a change of variable from $t$ to the true anomaly, using $dt=d\psi/\dot{\psi}$. The new limits of integration become $\psi_{i-1}$ and $\psi_{i}=\psi_{i-1}+2\pi$, or more simply $[0,2\pi].$ The orbital elements now depend on the true anomaly rather than time, which simplifies the integrands, but this is still not enough to evaluate them analytically. To do so, we use the fact that the orbits are osculating and the GW emission occurs instantaneously at pericenter, which ensures that $r_{p}$ and $e$ are constant everywhere except at closest approach. With this, the integrals become

[TABLE]

which can be evaluated analytically.

Alternatively, we can exploit the definition of orbital averaged GW fluxes to rewrite these changes in a simpler way. The orbital averaged energy flux, for example, is given by $\langle\dot{E}\rangle\equiv\Delta E_{(i,i-1)}/P$, where $P$ is the orbital period of the binary, and likewise for the angular momentum flux. With this definition, we are free to write

[TABLE]

Indeed, we recognize the integral expressions in Eqs. (14) and (15) as simply the product of the orbital period and $\langle\dot{E}\rangle_{i-1}$ or $\langle\dot{L}\rangle_{i-1}$ by definition. It might seem odd that orbit averaged quantities appear in the above expressions since the GW emission is happening mostly during pericenter passage, and thus smearing the emission over the entire orbit would appear incorrect. However, this is purely a result of the definition of the orbital averaged quantities, and has nothing to do with the nature of the GW emission or the validity of the orbital averaged approximation for the systems we are considering [60].

The orbital energy and angular momentum have a clear physical meaning, but $r_{p}$ and $e$ allow us to straightforwardly visualize the geometry of the system that is generating the bursts (at least at Newtonian order). At this order, it does not matter which set of quantities, $(E,L)$ or $(r_{p},e)$, we decide to use for the orbital evolution. For the Newtonian model in [19], we decided to use $(r_{p},e)$, so let us continue to do so here. We need to solve the system given by Eqs. (8) and (9) for the functionals $r_{p}\left(E,L\right)$ and $e\left(E,L\right)$. To obtain the evolution of the pericenter distance and the orbital eccentricity, we evaluate the functionals at the desired orbit, specifically $r_{p,i}=r_{p}\left[E_{i}\left(E_{i-1},L_{i-1}\right),L_{i}\left(E_{i-1},L_{i-1}\right)\right]$ and $e_{i}=e\left[E_{i}\left(E_{i-1},L_{i-1}\right),L_{i}\left(E_{i-1},L_{i-1}\right)\right]$. Evaluating the functionals with Eqs. (16) and (17) gives

[TABLE]

where we have kept only leading-order terms in $\delta e\equiv 1-e\ll 1$ and in $M/r_{p}\ll 1$ in all expressions (since we are working to Newtonian order). These equations recursively describe how the orbit shrinks and circularizes as the binary inspirals, thus completely describing the orbital evolution.

2.1.2 Centroid Mapping

The second ingredient we need for any burst model is the centroid mapping. The centroids of the bursts are given by the set $(t_{i},f_{i})$ at which the bursts occur. What we desire is the mapping $t_{i-1}\rightarrow t_{i}$ and $f_{i-1}\rightarrow f_{i}$. Since the orbits are osculating, the time between bursts is trivially given by the orbital period $P$ to Newtonian order, which is given by

[TABLE]

Thus to obtain the time mapping, we simply have to evaluate the orbital period at the desired orbit,

[TABLE]

where the mappings $r_{p,i}\left(r_{p,i-1},\delta e_{i-1}\right)$ and $\delta e_{i}\left(r_{p,i-1},\delta e_{i-1}\right)$ are given by Eqs. (18) and (19), respectively.

The fact that the bursts are separated by an orbital period can be seen more generally by writing $\dot{\psi}=\dot{\psi}_{\text{cons}}+\dot{\psi}_{\text{diss}}$, where $\dot{\psi}_{\text{cons}}$ is the conservative part coming from Keplerian orbital dynamics, and $\dot{\psi}_{\text{diss}}$ is the dissipative piece that comes from radiation reaction. The time between successive pericenter passages, i.e. the orbital period, is then

[TABLE]

The dissipative part contains terms that depend on the radiation reaction gauge, which vanish upon integration. Thus we are left with only the conservative piece of $\dot{\psi}$, and by assuming the orbits are osculating and the GW emission is instantaneous, this evaluates upon integration to the orbital period for an unperturbed, purely conservative orbit. This does not mean that the orbital period is not evolving. GW emission carries energy and angular momentum away from the binary that changes the orbital period. The dissipative part of the above integral does vanish, but the dissipative part of $P$, namely $\dot{P}$, does not. The above result simply implies that the time between pericenter passages is the orbital period of a Keplerian orbit, since radiation reaction is happening rapidly around pericenter.

The GW frequency on the other hand requires knowledge of the Fourier transform of the GWs emitted during each burst, the Fourier-domain waveform. From Fig. 7 in [61], the GW power is highly peaked around

[TABLE]

where $\tau_{\rm GW}$ is the characteristic GW time, defined by

[TABLE]

At Newtonian order, the pericenter velocity is given by

[TABLE]

and thus $\tau_{\rm GW}$ is

[TABLE]

which roughly corresponds to the amount of time the system spends at pericenter. This time is a functional of the pericenter distance and eccentricity; hence, to obtain the frequency of the $i$-th burst $f_{i}$, one simply has to evaluate the characteristic GW time associated with the orbit $(r_{p,i},e_{i})$, where the mapping to the parameters of the previous orbit are given by Eqs. (18) and (19):

[TABLE]

This completes the mapping of the time-frequency centroid of the bursts.

The prescription we provide for the frequency of the bursts is dependent on the frequency domain waveform, or alternatively the GW power, peaking at $\tau_{\rm GW}^{-1}$. This intuition comes from [61], where for parabolic orbits, and at Newtonian order, it is shown that the GW power peaks roughly at $\tau_{\rm GW}^{-1}$. However, for circular binaries, the power peaks at twice the orbital frequency, and it can easily be checked that the prescription given above does not reproduce this result when $e=0$. This implies that there are uncontrolled remainders that depend on $\delta e$ that correct the above expression to account for this. However, because we are working in the limit where $\delta e\ll 1$, we expect such corrections to be subdominant.

2.1.3 Volume Mapping

The last ingredient we need is the volume mapping. The bursts are not instantaneously emitted at pericenter and are not solely peaked at one frequency. The emission is instead spread out over the full pericenter passage and over multiple frequencies. To complete the burst model, we need to determine how the time-frequency size of the bursts change from one to another. We may describe the bursts as any two dimensional objects in time and frequency. For simplicity, we chose to model the bursts as boxes with widths

[TABLE]

where $\xi_{t}$ and $\xi_{f}$ are constants of proportionality that are chosen from data analysis considerations. For example, one can choose these constants such that a desired percentage of the GW power ($90\%$ for example) is contained in each box. More general two dimensional objects, such as ellipsoids, could be used for this construction, but boxes are the simplest. To obtain how these widths change from one burst to the next, we simply have to evaluate Eqs. (28)-(29) at the parameters of the orbit $(r_{p,i},e_{i})$:

[TABLE]

Note that in the case of ellipsoids, the results are the same, but these quantities can instead be interpreted as the semi-minor and semi-major axes of the ellipsoids. For a realistic search, ellipsoids would actually be more appropriate choice since they are a more accurate representation of the time-frequency structure of the bursts. However, for the purposes of this work, this choice is irrelevant, as the goal is to characterize the two dimensional objects via the scales in Eqs. (30) and (31). This completes the review of the construction of the burst model to Newtonian order.

2.2 A Simplified Formalism

Ultimately, we are interested in a PN burst model at generic (presumably very high) PN order. Building a generic order PN model by following the construction above might at first seem like an intractable problem. To start, one would have to take the orbital energy $E_{\rm PN}$ and angular momentum $L_{\rm PN}$ at arbitrary order and invert these expressions to obtain $r_{p}(E_{\rm PN},L_{\rm PN})$ and $e(E_{\rm PN},L_{\rm PN})$. Then, one would need to use the energy and angular momentum fluxes to compute the evolution of the orbital energy and angular momentum. From there, the pericenter and orbital eccentricity mapping would have to be computed using the functionals $r_{p,i}\left[E_{i}(E_{i-1},L_{i-1}),L_{i}(E_{i-1},L_{i-1})\right]$ and $e_{i}\left[E_{i}(E_{i-1},L_{i-1}),L_{i}(E_{i-1},L_{i-1})\right]$. While this may actually be possible from a mathematical standpoint, it will be very non-trivial to do so at arbitrary order. Thus, in this subsection, we will instead seek a simplified formalism that is more practical to implement.

The new method we seek must be more direct than the previous method discussed, removing steps that are redundant and reducing the number of physical quantities we need to work with. We begin by noting a number of assumptions that we will use to simplify the analysis:

1. I.

Osculating Orbits: Any changes in the orbital parameters will be modeled as occurring instantaneously around pericenter passage, leaving the orbital parameters constant throughout the rest of the orbit. 2. II.

High Ellipticity: The orbits we consider are highly elliptical, so we define a small parameter $\delta e_{t}\equiv 1-e_{t}$ and work perturbatively in the regime $\delta e_{t}\ll 1$. 3. III.

PN Orbits: We will work within the PN framework, expanding all expressions in the pericenter velocity $v_{p}\ll 1$.

The first and second assumptions follow directly from the nature of the systems we consider in this paper. Note that in the second assumption we are now working with the time eccentricity from the QK parametrization. The reason for this is that in PN theory, there is no unique concept for the orbital eccentricity, as there are actually three eccentricities that enter the QK equations of motion, specifically $(e_{t},e_{r},e_{\phi})$. All three of these eccentricities reduce to the orbital eccentricity in the Newtonian limit, but they are distinct quantities within PN theory. We choose to work with $e_{t}$ and will express all quantities in terms of it. We discuss this in more detail in B.

The final assumption is new to this analysis and replaces the previous Newtonian assumption. At Newtonian order, we worked with the pericenter distance $r_{p}$ as one of our physical parameters. We will now choose to work with the pericenter velocity $v_{p}$ instead. This change is meant to put the computation more in line with the standard PN formalism for quasi-circular inspirals, as well as to remove some difficulties that result in there being terms that depend on half-integer powers of $r_{p}$ in the dissipative sector. The mapping between the pericenter distance and velocity is given explicitly to 1PN order in Eqs. (163) and (165), with the 2PN and 3PN corrections given in Eqs. (B.2) and (B.2). With the above assumptions, any burst model requires the three ingredients laid out in Sec. 2.1.

Let’s start with the orbital evolution, where now we focus on the evolution of the pericenter velocity and the orbital eccentricity. Rather than starting from the orbital energy and angular momentum, we are free to write the velocity and eccentricity mappings as

[TABLE]

where $\Delta v_{p,(i,i-1)}$ and $\Delta e_{t,(i,i-1)}$ are the change in pericenter velocity and time eccentricity between two successive orbits, and we have used the fact that $\delta e_{t}=1-e_{t}$ to write $\Delta\delta e_{t}=-\Delta e_{t}$. These mappings are directly analogous to the mappings of energy and angular momentum in our Newtonian model, given by Eqs. (12) and (13).

Expressions for $\Delta v_{p,(i,i-1)}$ and $\Delta e_{t,(i,i-1)}$ can be found in exactly the same way as in Eqs. (14) and (15). We may thus jump ahead and directly write

[TABLE]

where $\dot{v}_{p}$ and $\dot{e}$ are the rates of change of pericenter velocity and orbital eccentricity. These rates, once again, depend on the true anomaly and thus have gauge-dependent terms arising from the GW sector. Upon integration, these terms vanish, except now the above quantities are not necessarily gauge-invariant as they depend on the specific coordinate system one chooses to do the PN calculation in.

To our knowledge, the expressions $\dot{v}_{p}\left(v_{p},e_{t},\psi\right)$ and $\dot{e}_{t}\left(v_{p},e_{t},\psi\right)$ have not yet been explicitly computed and would not be easy to compute, which would leave something of a gap in constructing the orbit evolution for our bursts. However, we may once again exploit the definition of orbit averaging and write Eqs. (32) and (33) as

[TABLE]

The orbit averaged quantities $\langle\dot{v}_{p}\rangle$ and $\langle\dot{e}_{t}\rangle$ can be easily computed from the orbital energy and angular momentum and the corresponding fluxes, which are known to full 3PN order. We have thus completely constructed the orbit evolution for our burst model.

We now focus on the centroid and volume mappings. Once again, we will treat the bursts as boxes in time and frequency, and determine the mapping between the centroids and widths of the boxes. The characteristic GW time is still given by Eq. (24). In the Newtonian model, we used the expression $v_{p}\left(r_{p},e_{t}\right)$ given by Eq. (25) to write this time in terms of $(r_{p},e_{t})$. Since we are now working with $v_{p}$ instead of $r_{p}$, we can invert the relationship between these two parameters to obtain $r_{p}\left(v_{p},e_{t}\right)$,which is given explicitly in Eqs. (163)-(165) and (B.2)-(B.2), and write $\tau_{\rm GW}$ in terms of $(v_{p},e_{t})$. Once this time is specified, we may define the characteristic GW frequency by Eq. (23). The centroid and volume mappings follow the exact same analysis as the Newtonian model, only parameterized by the pericenter velocity rather than the pericenter distance. Hence, we may write

[TABLE]

thus completing the last two ingredients we need for our simplified formalism.

3 A Generic PN Formalism

With the application of assumption I, we have constructed a purely generic burst model through Eqs. (36), (37), and (38)-(41) that applies in any theory of gravity. We now seek to use this formalism to create a burst model at generic PN order. We will provide explicit expressions for the burst model at 1PN, 2PN, and 3PN orders in Sec. 3.1.

The above considerations imply that, to construct our burst model, we need PN expansions for four quantities: the orbital period, the pericenter distance, the rate of change of pericenter velocity, and the rate of change of orbital eccentricity. We can write these expansions to arbitrary PN order as

[TABLE]

with the Newtonian order quantities

[TABLE]

We refer to the vector $\vec{X}$ as the PN state vector, which depends on the PN expansion parameter. In our case, the PN expansion parameter is $v_{p}$ and the components of $\vec{X}$ are simply $X_{k}=v_{p}^{k}$. Furthermore, we refer to the vector fields $(\vec{P},\vec{R},\vec{V},\vec{E})$ as PN amplitude vectors, which are functions of the orbital eccentricity and the symmetric mass ratio. In Eq. (42)-(45), we have chosen to include a $v_{p}$ label in the amplitude vectors to remind us that the functional form of its components depends on the parameter one expands about, i.e. if we had chosen to work with the $x$ PN expansion parameter instead of $v_{p}$, then the eccentricity and symmetric mass ratio dependence of the PN amplitude vectors would be different. The dot products between the state vectors and the amplitude vectors take the simple form

[TABLE]

where $\vec{A}\in(\vec{P},\vec{R},\vec{V},\vec{E})$. We recognize the above expression as the summation of PN corrections to the associated quantity. The summation index $k$ acts as the PN order of each term and starts at $k=2$ corresponding to the corrections at 1PN order. The components of $\vec{A}$, specifically $A_{k}(e_{t},\eta)$, we then recognize as the coefficient of the $k/2$-PN order term.

The components of the amplitude vectors $(\vec{P},\vec{R},\vec{V},\vec{E})$ can be easily computed from known PN quantities. As an example, consider the orbital period. This quantity can be written as a function of the reduced energy $\varepsilon$ and angular momentum $j$ through the equation $P=2\pi/n$, where $n$ is the mean motion, given to 3PN order by Eq. (348a) in [46]. In turn, the reduced energy and angular momentum can also be written in terms of the pericenter velocity and eccentricity, $\varepsilon(v_{p},e_{t})$ and $j(v_{p},e_{t})$, which can be inserted into the expression $P(\varepsilon,j)$ and expanded about $v_{p}\ll 1$. The coefficients of each power of $v_{p}$ are then the components of the vector field $P_{k}$. We will provide expressions for these components at specific PN orders when we construct burst models at specific PN orders.

This should be very reminiscent of the computations commonly carried out in the ppE formalism. In the latter, four deformations characterized by eight parameters (4 amplitude parameters, which are actually fields since they depend on the eccentricity of the binary, and 4 exponent parameters) are used to completely specify the burst model in a general theory of gravity that reduces to GR in the low-velocity/weak-field limit. Instead of the four amplitude parameters of the ppE model, we now have four amplitude vectors fields $\lambda_{\rm PN}^{a}=(\vec{P},\vec{R},\vec{V},\vec{E})$, which characterize the PN corrections to the Newtonian quantities. Also, the ppE exponent parameters have been replaced by the PN exponent $k$, which is a known number. Hence, instead of the eight ppE parameters, we need $4k$ PN functions when working to $k/2$-PN order. Each of the PN vector fields is a function of the parameters of the system, which we have written solely as functions of the eccentricity and the symmetric mass ratio. This will be true at 1PN order, but at higher PN order, the functions can depend on other physical parameters, such as the spins of compact objects, or the equation of state of supranuclear matter when at least one of the binary components is a NS.

The goal of this section will be to write the PN modifications to the Newtonian mappings in terms of the set of PN functions $\lambda_{{\rm PN},k}^{a}=(P_{k},V_{k},R_{k},E_{k})$. We begin with the first ingredient, the orbital evolution, specified in our simplified formalism by Eqs. (36) and (37). In particular, we concentrate first on the evolution of the pericenter velocity. By exploiting the definition of orbit averaging, we are able to write the change in this quantity as $P\;\langle\dot{v}_{p}\rangle$, which is exactly the second term in Eq. (36). Hence, to obtain the velocity mapping, we simply have to multiply Eqs. (42) and (44) together and expand in $v_{p}$. It is not difficult to see that our expansion is a product of two sums that is equivalent to a double sum of the form

[TABLE]

where we have used the definition of the Cauchy product to rewrite the product of the sums as the discrete convolution of two series. When $k-2<2$, the convolution is exactly zero. Using this result, we write the velocity mapping as

[TABLE]

where the Newtonian term in this expression is

[TABLE]

and the new amplitude vector $\vec{\cal{V}}$ is

[TABLE]

which should be expanded about $\delta e_{t}\ll 1$ by Assumption II. In the above expression, the pericenter velocity is decreasing from one orbit to the next for highly elliptic orbits at Newtonian order. We refer to this behavior as pericenter braking, which will be explored in more detail in Section 4.2.

We may follow the same procedure for the eccentricity mapping to find

[TABLE]

with the Newtonian function

[TABLE]

and the amplitude vector

[TABLE]

We thus find that the PN amplitude vectors $(\vec{\cal{V}},\vec{\cal{D}})$ can be expressed in terms of the known PN amplitude vectors $(\vec{P},\vec{V},\vec{E})$. The above expressions are purely generic within the PN formalism, allowing them to be applied at any PN order.

Now, let us consider the second ingredient of the PN burst model: the centroid mapping. The GW time is given in Eq. (24), while the pericenter distance is given in Eq. (43). We thus have that the GW time at arbitrary PN order is

[TABLE]

and the frequency mapping between boxes is

[TABLE]

where the functionals $v_{p,i}(v_{p,i-1},\delta e_{t,i-1})$ and $\delta e_{t,i}(v_{p,i-1},\delta e_{t,i-1})$ are given in Eqs. (53) and (57), respectively. The components of the amplitude vectors $\vec{R}^{(-1)}$ are defined recursively in A. The time mapping can trivially be constructed from Eq. (42) via

[TABLE]

which completes our calculation of the centroid mapping.

We have here not inserted Eqs. (53) and (57) into Eqs. (3) and (61), and re-expanded about the pericenter velocity being small to determine the time and frequency of the bursts. The reason for this is that such an expansion results in a significant loss of accuracy compared to numerical evolutions. This results from the behavior of the orbital period in the burst model which behaves as

[TABLE]

where $A$ is a constant. Since both $\delta e_{t,i-1}$ and $v_{p,i-1}$ are assumed to be simultaneously but independently small, expanding such a function about only one of them would impose an assumption on their ratio that is not justified.

Finally, we consider the volume mapping of the bursts. Once again, we treat the bursts as boxes in time and frequency with widths defined by Eqs. (28) and (29). Hence we simply have to evaluate these expressions within our PN formalism at $(v_{p,i},e_{t,i})$, thus obtaining

[TABLE]

[TABLE]

Not surprisingly, we see that one set of corrections, specifically the PN corrections to the pericenter distance, are the PN corrections to the frequency and box size mappings. This is exactly like in the ppE burst formalism, where one ppE deformation (with two parameters: $\beta_{\rm ppE},\bar{b}_{\rm ppE}$) characterized these mappings.

3.1 Example PN Burst Models

We have applied the fully general formalism of Section 2.2 to PN theory, developing a burst model at generic PN order. This arbitrary order model is characterized by four amplitude vector fields $\lambda^{a}_{\rm PNburst}=\left(\vec{P},\vec{R},\vec{\cal{V}},\vec{\cal{D}}\right)$, which can easily be constructed from the PN corrections to the orbital period, pericenter velocity, and rate of change of pericenter velocity and orbital eccentricity. These amplitude vectors are dependent on the set of GR parameters that characterize the system, namely $\lambda^{a}_{{\mbox{\tiny GR}}}=\left(\delta e_{t},\eta,...\right)$. We will now apply this formalism to generate a few example burst models at specific PN orders. There are multiple coordinate systems used to calculate PN quantities. Two that are typically used within the literature are the ADM and modified harmonic coordinates. We will choose to work within the ADM coordinates. The expressions in modified harmonic coordinates can easily be obtained through the appropriate coordinate transformations, which are given for example in Eq. (7.11) in [54].

3.1.1 Burst Model at 1PN Order

We begin by calculating the burst model to 1PN order. Recall that the model has three ingredients: the orbit evolution, the centroid mapping, and the volume mapping. We begin with the orbit evolution, which in our generic order model is given by Eqs. (53) and (57). There are no 0.5PN order corrections to any of the quantities considered here, so the state vector has only one component, specifically

[TABLE]

To achieve a burst model at 1PN order, we simply have to compute the 1PN functions $({\cal{V}}_{2},{\cal{D}}_{2})$. The functions $({\cal{V}}_{k},{\cal{D}}_{k})$ are given in general by Eq. (56) and (59), respectively. Setting $k=2$, these functions become

[TABLE]

where $(V_{2},E_{2},P_{2})$ are given in B. Working to ${\cal{O}}(\delta e_{t})$, the orbit evolution becomes

[TABLE]

with

[TABLE]

Let us now consider the centroid mapping. The evolution of the time centroid of the bursts is trivially given by the orbital period, so to 1PN order

[TABLE]

with

[TABLE]

where $v_{p,i}(v_{p,i-1},\delta e_{t,i-1})$ and $\delta e_{t,i}(v_{p,i-1},\delta e_{t,i-1})$ are given by Eqs. (69) and (70), respectively.

We now move onto the frequency centroid mapping, which is characterized by the functions $R_{k}^{(-1)}$. Using the recursion method is A, $R_{2}^{(-1)}=-R_{2}$, and the frequency centroid mapping becomes

[TABLE]

with

[TABLE]

Finally, we focus on the volume mapping, which is trivially given by the same corrections as the frequency centroid mapping:

[TABLE]

where we have defined

[TABLE]

This completes the burst model at 1PN order.

3.1.2 Burst Model at 2PN Order

Let us now calculate the burst model to 2PN order. The state vector has three components corresponding to 1PN, 1.5PN, and 2PN orders, specifically

[TABLE]

We begin by computing the orbital evolution in the burst model. To 2PN order, the pericenter velocity and eccentricity mappings become

[TABLE]

The Newtonian and 1PN order mappings do not change from the 1PN order model, and they are given in Eqs. (71)-(72) and Eqs. (69)-(70), respectively. Generally, the 1.5PN order and 2PN order components of the amplitude fields are given by

[TABLE]

Using the results of B, we obtain

[TABLE]

where we have used the results of [62] to evaluate the tail enhancement factors.

Next, let us consider the time centroid mapping, which at 2PN order is

[TABLE]

where the mappings $v_{p,i}(v_{p,i-1},\delta e_{t,i-1})$ and $\delta e_{t,i}(v_{p,i-1},\delta e_{t,i-1})$ are now given by Eqs. (3.1.2) and (3.1.2), respectively. There is no 1.5PN order correction to these expressions, since $P_{3}=0$. This is a result of the fact that the orbital period comes from the conservative orbital dynamics, and is thus symmetric under time reversal. Once again, the Newtonian time centroid mapping and 1PN order correction do not change from the 1PN order burst model. The 2PN order correction is characterized solely by $P_{4}(e_{t},\eta;v_{p})$, which is given in B. Expanding about $\delta e_{t}\ll 1$, we obtain

[TABLE]

Finally, consider the frequency centroid and box size mappings, which at 2PN order are

[TABLE]

The Newtonian and 1PN order terms are the same as those in the 1PN order burst model. Using the results of A, the field $R_{4}^{(-1)}$ is in general given by

[TABLE]

Applying the expressions for $(R_{2},R_{4})$ from B, we obtain

[TABLE]

This completes the burst model at 2PN order.

3.1.3 Burst Model at 3PN Order

Let us now extend the burst model to the current limits of our understanding of eccentric binaries within PN theory, i.e. to 3PN order. The state vector will now extend to $v_{p}^{6}$, specifically

[TABLE]

At 3PN order, the orbital evolution equations become

[TABLE]

The new functions $[{\cal{V}}_{5},{\cal{V}}_{6}]$ and $[{\cal{D}}_{5},{\cal{D}}_{6}]$ give the coefficients of the 2.5PN and 3PN order corrections of the orbital evolutions. In terms of the components of the amplitude vector fields $[\vec{V},\vec{E},\vec{P}]$, they are given by

[TABLE]

where we have used the fact that $P_{3}(e_{t},\eta;v_{p})=0=P_{5}(e_{t},\eta;v_{p})$. Using the results of B, we find for the 2.5PN order functions

[TABLE]

where we have used the results of [62] and neglected the 2.5PN memory terms. The constants $[\nu_{0},\nu_{1},\rho_{0},\rho_{1}]$ depend on the coefficients of the Padé approximants created for the 2.5PN order tail enhancement factors $[\psi(e_{t}),\tilde{\psi}(e_{t})]$ in [62]. The exact rational form of the coefficients are too lengthy to provide here. We simply give their numeric values, which are

[TABLE]

For the 3PN order functions, we find

[TABLE]

This completes the orbital evolution to 3PN order.

The time centroid mapping at 3PN order becomes

[TABLE]

where once again there is no 2.5PN order corrections since the orbital period comes from the conservative orbital dynamics. The 3PN order function $P_{6}(e_{t},\eta;v_{p})$ is given in B. Expanding about $\delta e_{t}\ll 1$, we obtain

[TABLE]

Finally, the frequency and box widths mappings at 3PN order are

[TABLE]

where the functions $R_{6}^{(-1)}(e_{t},\eta;v_{p})$ is

[TABLE]

Using the results in B, we obtain

[TABLE]

This completes the burst model at 3PN order.

4 Properties of the PN Burst Model

With the burst model complete to 3PN order, we complete this paper with some results that describe properties of the model. We begin by discussing the accuracy of the burst model when compared to numerical evolutions of the PN radiation reaction equations. Finally, we discuss a previously unreported phenomenon associated with the evolution of the pericenter velocity under radiation reaction.

4.1 Accuracy of the Burst Model

The burst model is meant to be an accurate representation of GW bursts emitted by highly elliptic binaries in nature. Further, since this model is designed to be used as a prior in data analysis for detecting such systems, it is paramount that we characterize the accuracy of the model. The ideal test of such an analytic model would be to compare the time of arrival and frequency of eccentric bursts from a numerical relativity simulation to the those from the burst model. However, there are currently no accurate numerical relativity waveforms for the highly elliptic systems considered here. Even the second best comparison, the same as above but with accurate PN waveforms, is also currently inapplicable due to the lack of such waveforms. With the two most ideal tests out of reach, we are left with comparing the burst model to the orbital evolution of binary systems (instead of their associated waveforms) under PN radiation reaction. Such a comparison allows us to gauge the accuracy of the approximations used to construct the burst model, as well as estimate the typical error we can expect when comparing to physically accurate waveform models.

We begin by describing the method through which we obtain the numerical evolution. Ideally, the equations we would want to numerically evolve are $\langle\dot{e}_{t}\rangle(v_{p},e_{t})$ and $\langle\dot{v}_{p}\rangle(v_{p},e_{t})$. However, as we will explain in Sec. 4.2, there is always a point $(v_{p},e_{t})$ where $\langle\dot{v}_{p}\rangle=0$ during the inspiral, which numerical routines will have difficulty integrating past. An alternative approach is to use a parameterization of the equations that does not present this behavior, e.g. $\langle\dot{e}_{t}\rangle(x,e_{t})$ and $\langle\dot{x}\rangle(x,e_{t})$. The expression for $\langle\dot{e}_{t}\rangle(x,e_{t})$ to 3PN order, neglecting memory contributions, is provided in Eqs. (6.18)-(6.19), (6.22), and (6.25) in [55]. To obtain the expression for $\langle\dot{x}\rangle(x,e_{t})$ to 3PN order, we follow the method detailed in B for $\langle\dot{v}_{p}\rangle$, which we summarize here. We begin by obtaining an expression for $x(\epsilon,j)$ by inverting Eq. (6.5) in [55]. We then take a time derivative and apply the chain rule, using the 3PN order expressions for the energy flux [54] and the angular momentum flux [55]. We expand the resulting expression in $x$ to obtain $\langle\dot{x}\rangle(x,e_{t})$.

For our numerical evolutions, we integrate the equations $\langle\dot{x}\rangle(x,e_{t})$ and $\langle\dot{e}_{t}\rangle(x,e_{t})$ including all of the instantaneous and tail contributions to 3PN order. For the tail enhancement factors, we use the analytic expressions provided in [62]. The initial conditions for the evolutions are set to guarantee the initial eccentricity is $e_{t,0}=0.9$ and the initial GW frequency is $f_{\rm GW,0}=10{\rm Hz}$, i.e. we use these initial conditions to solve for the initial value of $v_{p}$ using $f_{\rm GW}(v_{p},e_{t})$, which is provided in B. We then use the expression $v_{p}(x,et)$, which is obtained from the 3PN extension of Eq. (B.1) with Eq. (7.10) in [54], to obtain the initial value of $x$. For the three systems we study, a $(1.4,1.4)M_{\odot}$ NSNS binary, a $(1.4,10)M_{\odot}$ NSBH binary, and a $(10,10)M_{\odot}$ BHBH binary, the initial conditions are listed in Table 1. With the initial conditions set, we numerically integrate the equations using the NDSolve routine in Mathematica until we reach the time when

[TABLE]

which denotes the maximum value of $x$ for which test particle orbits are stable around a Schwarzschild black hole, i.e. we require that $p>2M(3+e_{t})$, where $p$ is the semi-latus rectum of the orbit and we have used the Newtonian relation $x=(M/p)(1-e_{t}^{2})$. Beyond this point, we consider the inspiral to be formally over and to use the burst model one would have to extend it to include merger and ringdown.

For the comparison to the burst model, we use $x(t)$ and $e_{t}(t)$ to construct the pericenter velocity as a function of time $v_{p}(t)$, which we then use with the results of B to obtain the orbital period and GW frequency as a function of time, specifically $P(t)$ and $f_{\rm GW}(t)$. To compute the values of these in the burst model to 3PN order, we start the model with the same initial conditions used for the numerical evolution. Once $(v_{p,0},e_{t,0})$ are specified, all future $(v_{p,i},e_{t,i})$ are determined from Eqs. (3.1.3)-(3.1.3). From here, the orbital period and the GW frequency are determined in the burst model from Eqs. (117) and (119).

Figure 1 shows the orbital period and GW frequency as functions of time in the burst model and the numerical evolution, as well as the relative error between the two. The relative error increases as time increases, but typically the error remains below $1\%$ for the first one hundred bursts. The reason the error increases is twofold. First, the eccentricity decreases as the binary inspirals due to the loss of energy and angular momentum by GW emission. The burst model uses an expansion about $\delta e_{t}\ll 1$, and it is thus most accurate in this regime. This error can be improved by going to higher order in $\delta e_{t}$ within the burst model if one wishes.

The second reason for the increasing error is that as the binary inspirals, the GW power becomes smeared over more of the orbit. As a result, the binary’s evolution resembles less a set of discrete steps. The burst model, which is only valid when $\delta e_{t}\ll 1$, hinges on the osculating behavior of highly eccentric orbits. This error is more difficult to control, but one way of improving it would be to match the evolution in the burst model to an evolution when the eccentricity is small. However, in this paper, we are only interested in highly elliptic orbits where this matching is unnecessary. Regardless, as the figure shows, the error between the burst model and the numerical evolution is sufficiently small that we can begin to test the burst model in idealized data analysis scenarios.

4.2 Pericenter Braking

Let us begin by recalling that within our generic PN burst model, the change of pericenter velocity to Newtonian order is given by

[TABLE]

where the function $V_{\rm N}(e_{t})$ is

[TABLE]

Notice from Eq. (54), which provides $V_{\rm N}(\delta e_{t})$, that to first order in $\delta e_{t}$, the above expression is negative and $v_{p}$ is thus decreasing. This seems counterintuitive considering what we know about quasi-circular binaries, i.e. as the orbital separation $r$ decreases, the orbital velocity $v$ increases, since $v$ and $r$ are inversely related by Kepler’s third law.

This behavior becomes more confusing when we consider the apocenter velocity $v_{a}$. Just as we can calculate $\langle\dot{v}_{p}\rangle$ using the method detailed in B, we may also compute $\langle\dot{v}_{a}\rangle$. Following this method, and working to Newtonian order, we have

[TABLE]

Notice that this expression depends on $V_{\rm N}(-e_{t})$, which is always positive. The apocenter velocity is thus always increasing as the binary inspirals, just as we would expect from quasi-circular binaries. As a result, the pericenter and apocenter velocities have very different behavior depending on the eccentricity of the system.

Let us try to understand this counter-intuitive behavior. The function $V_{\rm N}(e_{t})$ is a third order polynomial in eccentricity with an oscillating sign and with the coefficient of the ${\cal{O}}(e_{t})$ term greater than unity. This means that there will be a critical point $e_{t,{\rm crit}}$ where the function is zero, $\langle\dot{v}_{p}\rangle(e_{t,{\rm crit}})=0$, and due to the aforementioned behavior of the coefficients $e_{t,{\rm crit}}<1$. Let us solve for this critical point. To Newtonian order we find

[TABLE]

where we have defined

[TABLE]

The Newtonian expression for the critical eccentricity evaluates to $e_{\rm crit}^{\rm N}\approx 0.5557306$. Such a critical point also exists at 1PN order, except that now it is a function of the mass ratio and the pericenter velocity:

[TABLE]

where we have defined

[TABLE]

This function evaluates to $e_{t,\rm crit}^{\rm 1PN}\approx 0.5557306-(0.06536872\eta+0.3457145)v_{p}^{2}$. The overall effect of the 1PN term is to decrease the value of the Newtonian critical point, but there is no value of $v_{p}<1$ or $\eta\in(0,1/4)$ for which $e_{t,\rm crit}=0$ at 1PN order.

Why does this behavior occur physically? The answer to this question lies in circularization. As the binary inspirals, energy and angular momentum are radiated away in such a way that the orbital eccentricity decreases, making the binary more and more circular. For quasicircular binaries, $v_{a}=v_{p}+{\cal{O}}(e_{t})$, but for highly-elliptic binaries, $v_{p}\gg v_{a}$. As a highly-elliptical binary inspirals, $v_{p}$ and $v_{a}$ will approach the same value since the eccentricity approaches zero. However, if the eccentricity is above the critical value, the two velocities will not approach the same value if they are both increasing initially. Instead, circularization causes pericenter to brake when the eccentricity is above $e_{t,{\rm crit}}$, so that $v_{p}$ can approach $v_{a}$. In fact, one can easily show from Eqs. (126) and (128), that $v_{a}$ and $v_{p}$ obey the following conservation law at Newtonian order:

[TABLE]

To further display this behavior, we plot the streamlines of $(\langle\dot{v}_{p}\rangle,\langle\dot{e}_{t}\rangle)$333We have rescaled the values of $\langle\dot{v}_{p}\rangle$ and $\langle\dot{e}_{t}\rangle$ in these plots to exemplify the behavior of the streamlines. This does not changes the results of this section. at Newtonian order in Fig 2 and at 1PN order in Fig. 3. Notice that in both plots, systems with values of $(v_{p},e_{t})$ above the critical value of the eccentricity, which is represented by the dashed line in Fig. 2 and the dotted line in Fig. 3, display the pericenter braking behavior that appears in the burst model. On the other hand, the pericenter velocity for systems below the critical eccentricity is always increasing.

This pericenter braking behavior is not a property of the burst model per se, but rather it is inherited from the PN radiation-reaction equations. One may worry that this pericenter braking behavior may disappear if treating the problem exactly (for example, through a numerical treatment). The right panel of Fig. 3, however, shows a zoom of the streamlines at small velocities, where we see that the braking behavior persists. We thus conclude that it is unlikely that pericenter braking is a artifact of the PN expansion.

Figure (3) also allows us to compare the critical eccentricity computed at Newtonian order, at 1PN order and numerically. The latter is obtained by solving the 1PN expression for $\langle\dot{v}_{p}(e_{t})\rangle=0$ to find $e_{t,\rm crit}$. As expected, the numerical inversion disagrees with the its 1PN expansion at high velocities. We notice, however, that the 1PN expression is closer to the numerical inversion than the Newtonian expression is. If the numerical inversion is correct, then this implies the 1PN expansion of $e_{t,\rm crit}$ given in Eqs. (131)-(4.2) has a larger regime of validity than its Newtonian counterpart.

Finally, it is important to note that while the pericenter velocity has this unique behavior, the GW frequency and the PN parameter $x$ are both monotonically increasing, and the time eccentricity is monotonically decreasing, throughout the inspiral of the binary. In the circular case, there is a one-to-one mapping between the orbital velocity and the GW frequency, and since the orbital velocity is a monotonic function, so is the frequency. For generic eccentric inspirals, the frequency depends on both the pericenter velocity (or alternatively $x$) and the time eccentricity in such a way that it is also monotonic.

5 Discussion

We have constructed a generic PN order burst model. This model is characterized by four amplitude vector fields $(\vec{P},\vec{R},\vec{\cal{V}},\vec{\cal{D}})$, which depend on the orbital period, pericenter distance, and rates of change of pericenter velocity and orbital eccentricity, respectively. While these quantities are not typically reported within the literature, they can be easily calculated from the quantities that are. Thus, the formalism presented here provides a formulaic means of generating burst models to any PN order. We have then applied this formalism to calculate the burst model out to the current limit to which we can compute PN quantities for eccentric binaries, i.e. 3PN order.

One direction of future research is to relax some of the assumptions used to develop this formalism. For example, we have approximated the compact objects as non-spinning point particles, which is appropriate if we are considering non-spinning BHs. However, BHs in the universe are generally considered to be spinning, while on the other hand, NSs are not well approximated by point particles. NSs will typically have small spins, however the inclusion of finite size effects and tidal perturbations would be necessary to effectively model highly elliptic NS binaries. Further, if one of the binary components is a BH, then not all of the GW power travels to spatial infinity. Instead, some of the GWs travel through the horizon of the BH, increasing its mass and spin throughout the evolution of the binary. With these considerations, we can postulate that the generic PN formalism can be extended to include such effects by writing

[TABLE]

where $\vec{A}\in(\vec{P},\vec{R},\vec{\cal{V}},\vec{\cal{D}})$. In the above, $\vec{A}_{\rm PP}$ represents the point particle terms, computed here to 3PN order, $\vec{A}_{\rm Spin}$ are the corrections generated by the spins of the compact objects, $\vec{A}_{\rm FS}$ are generated by finite size effects of NSs, and $\vec{A}_{\rm H}$ incorporates the corrections from the GWs fluxes through BH horizons. The final term, $\vec{A}_{\rm ppE}$ represents corrections due to modified theories of gravity, which have already been considered in [19].

One important question to address in the future concerns the most appropriate equations one should use to obtain the numerical evolution of highly elliptic systems under radiation reaction. In this work, we have used the orbit averaged equations for $\langle\dot{e}_{t}\rangle$ and $\langle\dot{x}\rangle$. These equations are applicable when the GW emission is smeared over the entire orbit and changes to the orbital elements are small on the timescale of one orbit, as is the case in quasi-circular inspirals. However, for the highly elliptic binaries considered here, the GW emission is concentrated at pericenter passage, and changes to the orbital elements happen on timescales significantly shorter than the orbital period. The evolution of such binaries will resemble a set of discrete steps from one orbit to the next. Furthermore, it can be shown that when expressed in terms of variables that are finite in the parabolic limit, the orbit averaged fluxes of energy and angular momentum vanish for parabolic orbits. There is, of course, nothing special about the parabolic limit, and binaries on parabolic orbits will still emit GWs, which suggests a break down of the orbit-averaged formalism in this limit. Since the orbit averaged equations are currently used prolifically in the literature, it is important to determine how big of a deviation in observables is generated by considering evolutions with and without orbit-averaging in the radiation-reaction force, and what set of systems in $(f_{\rm GW},e_{t})$ space are affected by this deviation. Such a study is currently underway [60].

Another avenue for future research is to consider how the 3PN order burst model aids in detecting highly elliptic binaries. In such a study, one would inject a waveform generated by numerically evolving the binary under radiation reaction into a simulated LIGO data stream. One could then perform an analysis to study whether the prior, specifically the burst model, is sufficient to achieve detection of such a signal given a particular noise model. One could also investigate the nature of posterior probability densities of recovered parameters and determine if such a search is accurate enough to perform parameter estimation on actual signals. With such a study completed, a follow up study could be conducted to investigate the search strategy’s ability to estimate deviations from the current model, such as those from modified theories of gravity and to constrain the coupling constants of such theories. Such studies will be crucial for understanding our ability to detect and perform important astrophysics with eccentric GW signals.

Acknowledgements

We would like to thank Frans Pretorius for several useful discussions. N.Y. acknowledges support from the NSF CAREER Grant PHY-1250636. N. L. acknowledges support from the NSF EAPSI Fellowship Award No. 1614203.

Appendix A PN Recursion Relations

When computing the burst model to arbitrary PN order, we are often faced will expressions of the form

[TABLE]

which need to be perturbatively expanded about $x\ll 1$. In our PN burst model, $v_{p}$ takes the place of $x$. The expansion of the above expression can be easily computed term by term, but may not be expressible in terms of an arbitrary sum. We instead define the coefficients as

[TABLE]

where it is understood that we are working perturbatively in $x$. If one computes the Taylor expansion and calculates the coefficients, one finds

[TABLE]

where the coefficients $(B_{n},C_{n},D_{n})$ are given by

[TABLE]

We have stopped the expansion at fourth order in $x$, but in principle there will also be fifth order terms with coefficients $E_{n}$, sixth order terms with coefficients $F_{n}$, etc., where ($E_{n}$, $F_{n}$, …) are expressible in terms of the coefficients above. Thus, while the right-hand-side of Eq. (A) is not expressible as a closed sum, one can recursively build the coefficients to arbitrary order.

As an example of how this works, consider the case where $m=1$ and we truncate at third order, i.e. the series only has coefficients $(A_{1},A_{2},A_{3})$. So, the first order term ($n=1$) is then

[TABLE]

However, from Eq. (138)

[TABLE]

Likewise, all higher order terms, i.e. $(C_{1},D_{1},...)$, will vanish, and thus at first order $A_{1}^{(-1)}=-A_{1}$ as one would expect from a first order Taylor expansion. At second order $n=2$,

[TABLE]

However, now the $B$ coefficient doesn’t vanish, but the $C$ coefficient (and thus all higher order coefficients) does vanish.

[TABLE]

So, the second order coefficient is $A_{2}^{(-1)}=-A_{2}+A_{1}^{2}$. Similarly, a calculation of the third order coefficient gives $A_{3}^{(-1)}=-A_{3}+2A_{1}A_{2}-A_{1}^{3}$. Since we are working perturbatively to third order in this example, we now truncate the series and all higher order coefficients will be zero.

The expression in Eq. (A) can be written in a more compact form using the notation of Section 3. Using the vector convolution $\circ$, we may write

[TABLE]

Appendix B PN Vector Fields

We here give explicit expressions for the PN amplitude vector fields $(P_{k},R_{k},V_{k},E_{k})$ to 3PN order. We will provide an example calculation at 1PN order to show how to obtain these functions. At higher PN order, we will simply give the results of the calculation.

B.1 1PN Amplitude Vector Fields

We begin by calculating the 1PN corrections to the factors $(P_{2},R_{2},V_{2},E_{2})$. This naturally amounts to finding the 1PN corrections to the orbital period, pericenter distance, and rates of change of pericenter velocity and orbital eccentricity. Because these are not quantities that are typically computed in the PN literature, we will show explicitly how to calculate them from other known quantities.

Consider first the pericenter velocity and the orbital eccentricity as functions of the orbital energy and angular momentum. Typically, PN quantities are written in terms of the reduced orbital energy $\varepsilon=-2E/\mu$ and the dimensionless orbital angular momentum $j=-2EL^{2}/\mu^{3}M^{2}$, where $\mu$ is the reduced mass of the system. The reason for this is that $\varepsilon$ and $j$ do not depend on the coordinate system that one chooses to do calculations in, i.e. they are coordinate invariant quantities. However, we have formulated our burst model in terms of the pericenter velocity and orbital eccentricity, so we need to determine the mappings $(\varepsilon,j)\rightarrow(v_{p},e)$.

We start from the equations of motion that govern the QK representation, which to 1PN order are [46]

[TABLE]

The first of these is the radial equation of the elliptical orbit written in terms of the eccentric anomaly $u$, where $a_{r}$ is the semi-major axis of the ellipse and $e_{r}$ is the radial eccentricity. The second equation is Kepler’s equation, which relates the eccentric anomaly to the mean anomaly through the time eccentricity $e_{t}$. This equation is itself a direct measure of time since $\ell=n(t-t_{p})$ where $n=2\pi/P$ is the mean motion and $t_{p}$ is the time of pericenter passage. Finally, the last equation is the azimuthal equation of the orbit, which relates the orbital phase $\phi$ to the eccentric anomaly through the azimuthal eccentricity $e_{\phi}$ and the advance of periastron per orbit $K$. In Newtonian gravity, $K=1$ and we recover the Newtonian equations of motion, but at 1PN order $K=1+{\cal{O}}(\varepsilon)$. Thus a primary difference in the azimuthal motion at 1PN order is the inclusion of the precession of periastron.

Another difference between Newtonian and 1PN orbits in the QK parametrization is the need for multiple ”eccentricities.” The latter are not actually separate physical quantities, but instead are a clever notational trick that allows the equations of motion to take the same functional form as their Newtonian analogs444Note that while this is true to 1PN order, the same cannot be said at higher PN order. At 2PN order and beyond, Eqs (148) and (149) pick up higher harmonics of the anomalies.. In reality, these eccentricities are related to the Newtonian expression for the orbital eccentricity, which is just $e=\sqrt{1-j}$, through [46]

[TABLE]

In this paper, and in this Appendix, we choose to parameterize the motion in terms of one of these eccentricities, specifically $e_{t}$, and we define $\delta e_{t}=1-e_{t}$. This choice is arbitrary: one could easily construct the burst model in terms of $e_{r}$ or $e_{\phi}$ or the Newtonian $e$. With this choice, we now have one of the equations we need to determine $\varepsilon(v_{p},e_{t})$ and $j(v_{p},e_{t})$, specifically Eq. (151).

Let us now focus on the pericenter velocity. Since pericenter is the minimum turning point of the orbit ($\dot{r}=0$), we may write $v_{p}=r_{p}\dot{\phi}(r=r_{p})$. From Eq. (147), we see that $r$ is at a minimum when $u=0$, specifically $r(u=0)=r_{p}=a_{r}(1-e_{r})$. Hence we are left with determining $\dot{\phi}(u)$. Taking the time derivative of Eq. (149), we find

[TABLE]

where $\beta=\sqrt{(1+e_{\phi})/(1-e_{\phi})}$. To find the expression for $\dot{u}$, we take the time derivative of Eq. (148) and solve for $\dot{u}$,

[TABLE]

We can now put all of this together to find $v_{p}$ to 1PN order:

[TABLE]

This expression for the pericenter velocity is exact in the sense that we have not performed a PN expansion in $\varepsilon$ yet. To obtain $v_{p}(\varepsilon,j)$, we use Eqs. (150)-(152) with [46]

[TABLE]

and expand in $\varepsilon$ to find

[TABLE]

The last step before deriving the functions $(P_{2},R_{2})$ is to invert Eqs. (151) and (B.1) to find $\varepsilon(v_{p},e_{t})$ and $j(v_{p},e_{t})$. This can be done very easily order by order in $\varepsilon$ (or $v_{p}$) to find

[TABLE]

We can now construct the functions $P_{2}$ and $R_{2}$, which recall we define via

[TABLE]

The orbital period in the QK representation is given by $P=2\pi/n(\varepsilon,j)$, with $n(\varepsilon,j)$ given in Eq. (156). Inserting our expressions for $\varepsilon(v_{p},e_{t})$ and $j(v_{p},e_{t})$ and expanding about $v_{p}\ll 1$, while keeping terms of relative ${\cal{O}}(v_{p}^{2})$, we find

[TABLE]

We may follow the same procedure for the pericenter distance, using $r_{p}=a_{r}(\varepsilon,j)\left[1-e_{r}(\varepsilon,j)\right]$; we find

[TABLE]

The functions $(V_{2},E_{2})$ require a bit more work, since we need to compute the rates of change of pericenter velocity and time eccentricity. We may do this via the chain rule:

[TABLE]

where $h=L/M$ is the reduced angular momentum, $\dot{E}$ is the energy flux, and $\dot{L}$ is the angular momentum flux. The reason for replacing $j$ with $h$ in the expressions $\left[v_{p}(\varepsilon,j),e_{t}(\varepsilon,j)\right]$ is that $j$ has a factor of the reduced energy hidden in it, which needlessly obfuscates taking partial derivatives. Furthermore, the expression for $\dot{h}$ can be computed directly and easily from $\dot{L}$. The averaged energy and angular momentum fluxes are given to 1PN order by [54, 55]

[TABLE]

with the enhancement factors

[TABLE]

and

[TABLE]

In the above, the PN expansion parameter $x=(M\Omega)^{2/3}$, with $\Omega$ the orbital frequency, which is related to the reduced orbital energy by [46]

[TABLE]

The goal at this stage should be clear: we desire to write $\langle\dot{v}_{p}\rangle$ and $\langle\dot{e}_{t}\rangle$ in terms of $v_{p}$ and $e_{t}$ instead of $x$ and $e_{t}$. The calculation is rather lengthy, but straightforward, so we will simply detail the main steps here. The first step is to replace $x$ with $\varepsilon$ in the expressions for $\langle\dot{E}\rangle$ and $\langle\dot{L}\rangle$, working perturbatively in $\varepsilon$. We then use the chain rule combined with our new expressions for the fluxes in Eqs. (166) and (167). Finally we replace any instances of $(\varepsilon,h)$ with $(v_{p},e_{t})$ and expand in powers of $v_{p}$. The end result is

[TABLE]

where $\langle\dot{v}_{p}^{\rm N}\rangle$ and $\langle\dot{e}^{\rm N}\rangle$ are given in Eqs. (48) and (49) with the replacement $e\rightarrow e_{t}$, and the PN functions are

[TABLE]

This completes the construction of the 1PN correction functions $(P_{2},R_{2},V_{2},E_{2})$.

B.2 Amplitude Vector Fields to 3PN Order

We now provide the components of the amplitude vector fields $(\vec{P},\vec{R},\vec{V},\vec{E})$ to 3PN order. The derivation of the components follows the exact same procedure as the previous section, so we will simply list the components, where recall we work in ADM coordinates. For the orbital period, the non-zero components are

[TABLE]

For the pericenter distance, the non-zero components are

[TABLE]

For $\langle\dot{e}_{t}\rangle$, we have

[TABLE]

[TABLE]

where $\gamma_{E}=0.577...$ is the Euler constant and the tail enhancement factors $[\varphi_{e}(e_{t}),\psi_{e}(e_{t}),\zeta_{e}(e_{t}),\kappa_{e}(e_{t})]$ are defined in [55]. The functions $[{\cal{P}}_{e}^{(k)}(\eta),{\cal{S}}_{e}^{(k)}(\eta)]$ are polynomials in $\eta$, specifically

[TABLE]

Finally, for $\langle\dot{v}_{p}\rangle$, we have

[TABLE]

where $[\phi(e_{t}),\psi(e_{t}),\zeta(e_{t}),\kappa(e_{t})]$ and $[\tilde{\varphi}(e_{t}),\tilde{\psi}(e_{t}),\tilde{\zeta}(e_{t}),\tilde{\kappa}(e_{t})]$ are the tail enhancement factors for the energy and angular momentum fluxes, respectively. These functions are defined in  [63, 55], with analytic representations as functions of $e_{t}$ given in [62]. The polynomials $[{\cal{P}}_{v}^{(k)}(\eta),{\cal{S}}_{v}^{(k)}(\eta),{\cal{L}}_{v}(e_{t})]$ are given as follows,

[TABLE]

This completes the amplitude vector fields to 3PN order.