is the overall waveform (including upsampling). With similar requirements to
EMRIs given in the section above, we expect to collect at least $\sim 10^{6}$
MCMC samples for a given source posterior distribution. To accomplish this in
$\sim 10$ days, we need to generate waveforms at a rate of $\sim/1$s on a
single central processing unit (CPU) core.

The other efficiency requirement deals with the separation of the accurate
waveform generation from the scaling operation. Generally, the scaling
operation is the bottleneck for LISA where waveforms can have up to $\sim
10^{6}$–$10^{8}$ data points. Therefore, the requirement on the scaling part
is that a waveform on a single CPU core be scaled in $\sim 1$ second or less.
This condition, therefore, also sets the requirement for the accurate waveform
portion assuming its sparsity prevents it from becoming the bottleneck. This
requirement is that the sparsely sampled, accurate waveform must be of a
similar or lower order of magnitude in timing when compared with the scaling
operation.

###  3.4 Interface and data format requirements

Contributors: Tyson Littenberg

An additional key consideration for optimally supporting LISA analyses is the
need for a flexible interface between waveform generation software and
analysis pipelines. It is important for verification and validation of
pipelines, and as a means for cross-checking results, to have independently
developed algorithms targeting the same sources. (There is also the added
benefit of a constructive competition between development teams.) As per the
previous discussion around Eq. (4), different analysis pipelines will likely
be built around different representations of the data. Template generation is
generally the computational bottleneck, and so it needs to be optimized for
the application. At the same time, the benefits of having independent
pipelines become liabilities if the analyses are not interfacing to the same
waveform generation tools. As a result, it is of paramount importance that the
waveform and pipeline development teams are collaborating early and often to
avoid unnecessary or redundant transformations of the waveforms by the
analysis pipelines (consider, e.g., a waveform that is initially computed in
the time domain but then output in the Fourier domain, being called by an
analysis pipeline that uses a discrete wavelet domain representation of the
data).

The demand for flexibility of the waveform-analysis interface affects more
than just the choice of a basis set used to represent the template. There have
been promising developments in low-cost likelihood evaluations that use the
instantaneous amplitude and phase of the template waveform, sampled on an
adaptive grid, to concentrate computations to regions where the signal is
changing most rapidly [502, 504]. Such considerations are difficult to
retroactively incorporate in established template generation algorithms, but
present opportunities for increased efficiency if they are part of the
original waveform algorithm designs.

##  4 Modelling approaches for compact binaries

Detecting and inferring the parameters of compact binaries systems requires
waveform template which, in turn, necessitates solving the relativistic two-
body problem. Unlike in the Newtonian counterpart, it is not possible to only
solve for the motion of the two bodies; one must solve for the dynamical
evolution of the full spacetime. In general there are no known closed-form
solutions to the nonlinear Einstein field equations for radiating binaries and
so a variety of techniques have been developed to compute solutions either
numerically or via perturbative expansions. The three main approaches that
directly solve the Einstein field equations are numerical relativity (NR),
post-Newtonian/Minkowskian (PN/PM) theory, and gravitational self-force (GSF).

Each of these approaches has strengths and weaknesses that leads them to being
best employed in different regions of the binary configuration parameter
space. Numerical Relativity directly solves the Einstein field equations to
produce exact solutions up to numerical error. Post-Newtonian theory
analytically computes relativistic corrections to the binary’s motion and GW
emission as a series expansion in powers of the orbital velocity as a fraction
of the speed of light. The closely related post-Minkowskian approach expands
field equations around flat Minkowski spacetime in powers of the gravitational
constant without any restriction of the velocity of the binary. Gravitational
self-force expands the Einstein field equations in powers of the (small) mass
ratio. Fig. 2 gives a quantitive description of the strengths and weakness of
each approach in the orbital separation – mass-ratio parameter space for non-
spinning quasi-circular binaries.

In addition to the above approaches there are also effective frameworks that
attempt to cover large portions of the parameter space. The physically
motivated Effective-One-Body model takes inspiration from the solution to the
Newtonian problem and describes the binary as the motion of a test body in
spacetime of a deformed single black hole. Much progress can be made
analytically with this approach by absorbing post-Newtonian and post-
Minkowskian corrections, and further calibration can be applied using
numerical results from the NR and GSF methods. The Phenomenological (Phenom)
waveform models do not attempt to solve the relativistic two-body problem.
Instead they aim to directly model the waveforms by building upon fast post-
Newtonian models with further calibration from NR and GSF. The EOB and Phenom
models are heavily used in analysis of GW data from current ground-based
detectors.

In this section we outline the status of the above approaches and discuss the
required development needed to reach the accuracy requirements outlined in
Sec. 3.2. Concerns about the speed of waveform generation are addressed in
Sec. 5. The status and requires for LISA of these approaches can be found in
Sec. 4.1 for numerical relativity, Sec. 4.2 for weak-field post-Newtonian and
post-Minkowskian expansions, Sec. 4.3 GSF, Sec. 4.5 for EOB, and Sec. 4.6 for
Phenomenological models. These sections focus on modelling within GR; see Sec.
6 for a discussion on modelling in alternate theories of gravity.

Figure 2:  Region of applicability of different approximation techniques for
non-spinning quasi-circular binary BH inspiral. The 1PA region is a prediction
derived from fitting NR data to an expansion of the form (14). The shaded
regions indicate ranges within which the cumulative orbital phase-error is
less than $\pi/4$ and $\pi/16$ radians, respectively. Recent direct
calculations of post-adiabatic (1PA) waveforms suggest the 1PA region shown is
over-optimistic for $\nu>0.1$, but they have borne out the prediction for
$\nu<0.1$ [505]. The gray lines show the location of binaries with 10 (resp.
100 and 1000) orbits left before they reach the ISCO. This figure is
reproduced from Ref. [506].

###  4.1 Numerical relativity

Coordinators: Mekhi Dhesi, Deborah Ferguson, and Deirdre Shoemaker
Contributors: S. Bernuzzi, G. Bozzola, K. Clough, D. Ferguson, W. Gabella, M.
Gracia, R. Haas, M. Hannam, E. Huerta, S. Husa, L. Kidder, P. Laguna, C.
Lousto, G. Lovelace, D. Neilsen, V. Paschalidis, H. Pfeiffer, G. Pratten, H.
Rüter, M. Ruiz, S. Shapiro, J. Thompson, A. Tsokaros, H. Witek, M. Zilhao

####  4.1.1 Description

Numerical relativity solves Einstein’s field equations through direct
numerical integration on supercomputers (see e.g. the books [507, 508, 509,
510, 511]), providing the spacetime geometry and dynamics of the system in
addition to the gravitational radiation emitted. Such solutions of the full
nonlinear field equations without approximations are critical to our
understanding of highly dynamic regimes without symmetries, such as the late
inspiral and merger of compact binaries.

The first steps at numerical solutions of Einstein’s equations date back
several decades [512, 513] and tremendous progress since then enabled the
first simulations of inspiraling and merging binary black holes [514, 515,
516]. These breakthroughs initiated an explosive growth of the field, as
illustrated, e.g. in the reviews [517, 518, 519, 520]. By now, NR has become a
critical part of the modeling the gravitational waves of the late inspiral and
merger phase of coalescing binaries (e.g. [521, 522, 523, 524]), which in turn
underpins the analyses of all observed gravitational waves from coalescing
binaries.

In order to solve Einstein’s field equations, modern NR codes generally use
the 3+1 approach [525, 526] in which the four-dimensional spacetime is sliced
into three-dimensional hypersurfaces, and Einstein’s equations are
reformulated as a Cauchy problem with constraints. The first stage in
numerically solving a compact binary inspiral is the construction of initial
data, for which Einstein’s constraint equations are reformulated as elliptic
equations either in the context of the conformal-thin-sandwich formulation
[527, 528] or the puncture approach [529]. The resulting coupled nonlinear
partial differential equations are solved with custom-purpose elliptic solvers
[530, 531, 532]. Numerous improvements over the years increased the generality
of the physical conditions that can be achieved and the numerical quality of
the solution, e.g. [533, 534, 535, 536, 537]. Moreover, entirely new codes
were developed [538, 539, 540, 541, 542]. The time evolution is encoded in a
set of coupled, hyperbolic partial differential equations with suitable gauge
conditions and boundary conditions, where three main approaches have emerged:
the generalized-harmonic formulation [543, 544, 545], the Baumgarte-Shapiro-
Shibata-Nakamura formulation [546, 547] and different versions of the Z4
formulation of Einstein’s equations [548, 549, 550, 551, 552, 553, 554] (see
Table 7). Gravitational radiation must be propagated to future null infinity
(see review [555]), which is often accomplished by extrapolation of GW modes
extracted at finite radius [556, 557]. Alternatively, Cauchy-Characteristic
Extraction [558, 559, 560] directly yields the waveforms at infinity, with
improved recovery of non-oscillatory GW modes and GW memory effects [561].
Gauge conditions and transformations at future null infinity must be treated
with care [562] to yield well-behaved numerical waveforms with well-defined
waveform modes [563]. Cauchy codes have difficulty resolving the memory effect
[523, 561]. Instead, most approaches to computing the memory effect use Cauchy
NR waveforms (or NR calibrated waveform models) and the asymptotic Einstein
equations to determine the unresolved memory effect in the Cauchy simulations
that is required to satisfy Einstein’s equations [564, 565, 566, 567, 568,
569]. This approximate approach agrees with the Cauchy-Characteristic
extracted waveforms, though comparisons over a wider range of the BBH
parameter space would be useful.

Once BBH simulations became possible, the NR community quickly achieved many
firsts, including the first simulations of unequal-mass BBH coalescences [570,
571], the first with spinning binaries [572], the first eccentric mergers
[573] and first comparisons with PN calculations [574, 575, 576, 577]. NR
calculations of the merger revealed features of the non-linear regime, among
them that BHs with spins aligned with the orbital angular momentum merge at
higher frequency [572], as well as the calculation of the recoil velocity
imparted on the remnant BH [570, 578, 571]; in particular it was found that BH
spins oriented approximately parallel to the orbital plane can lead to recoil
velocities of several 1000’s km/s [579, 580, 581].

Analysis of GW observations requires waveforms that cover the entire frequency
band of the relevant detector at high accuracy. This has motivated large
efforts to improve the accuracy of NR simulations [577, 582, 545, 583], length
of the inspiral [521, 584], and coverage of increasingly large portions of
parameter space in increasing detail. Parameter space exploration was
initially performed through community-wide Numerical INJection Analysis
(NINJA) [585, 586] and Numerical-Relativity-Analytical-Relativity (NRAR) [587]
collaborations, which also yielded important cross-checks between different
codes [588, 587]. Newer and more extensive parameter space surveys are listed
in Table 7. Waveform models developed with NR information are described in
Sec. 4.5 and 4.6. NR surrogate models [589, 590] are directly built on NR
simulations, and are particularly important for the analysis of high-mass BBH
systems like GW150921 [114, 115]. NR simulations can also be directly used for
GW data-analysis, where they serve as synthetic signals [591] for quantifying
the response of GW search and parameter-estimation pipelines [592, 593], and
to conduct indirect analyses of observations [594, 595, 596].

NR also offers important information for the ringdown phase, characterized by
an exponential decay as the remnant BH settles into a Kerr black hole. NR
results determined the remnant parameters (mass, spin, recoil velocity) [597,
598, 599, 600, 601, 602, 590]. Detailed NR calculations of the emitted
gravitational waves during merger and ringdown yield the initial amplitudes
and phases of quasi-normal ringdown modes and underpin theoretical studies of
what information about nonlinear processes are accessible through the emitted
gravitational waves [603], and how the ringdown phase can be used to test GR
and to probe the no-hair and area theorems [604, 605, 606, 607, 608, 609, 610,
611].

Due to our focus on the role of NR in the LISA mission, this section primarily
covers vacuum spacetimes with a short discussion on environmental effects that
include matter; the important work of NR with neutron stars is not included.

####  4.1.2 Suitable for what sources?

The coalescence of BBHs is a primary source for LISA for a vast range of BH
mass. Those BBH systems that merge in the LISA band require NR to produce the
waveforms during the late inspiral and merger. In order to quantify the mass
ranges that merge in the LISA band, we compute the masses for which the
frequency of ringdown and inspiral are both in the LISA band. The dominant
quasi-normal mode for a non-spinning black hole has a frequency of
$f=12.07\mathrm{mHz}\left(M/10^{6}M_{\odot}\right)^{-1}$ [612], which will be
within the LISA band for BH masses of roughly $10^{4}-10^{8}M_{\odot}$.
Turning to the inspiral, NR simulations of BBH typically cover a few tens of
orbits before merger, starting at GW frequencies of $f\sim
1\mathrm{mHz}\left(M/10^{6}M_{\odot}\right)^{-1}$, placing this part of the
inspiral into LISA band for $10^{3}-10^{7}M_{\odot}$. This estimate implies
that BBHs of mass $10^{4}-10^{8}M_{\odot}$ will have at least part of its
waveform within coverage of NR. NR is capable of modeling these binaries for
comparable mass ratios of less than 1:100.

NR simulations become increasingly costly with increasing mass-ratio and with
very large spins, whereas there is only minor dependence of computational cost
on other BBH parameters like spin-direction and orbital eccentricity.
Parameter space coverage is improving over time. The next section describes
what parameters NR currently covers and in Sec. 4.1.5 we discuss where efforts
are needed to achieve complete coverage of the potential range of source
parameters.

####  4.1.3 Status of the Field

Since the breakthrough in NR in 2005/2006 [514, 515, 516] and the subsequent
“gold-rush” to explore nonlinear phenomena in black hole and neutron star
mergers, NR has now matured into a reliable tool for accurately computing
gravitational waveforms needed to characterize GW data. In this section we
briefly describe the landscape of current and future numerical relativity
codes, their cross-validation and the parameter space coverage of currently
available numerical waveforms.

The community has developed a number of successful NR codes that target a
variety of goals for BBH spacetimes including covering increasing fractions of
the parameter space (so far mostly comparable mass ratios and moderate spins),
increasing number of GW cycles, accuracy of waveforms and interactions with
matter. Table 7 provides a comprehensive list of currently available NR codes,
some of which are briefly described in A. Their capabilities have been
sufficient for the detection and characterization of GW signals with current
ground-based GW detectors. However, data collected with future space- or next
generation ground-based instruments will have a higher SNR and its
interpretation requires much improved NR waveforms.

Code | Open | Public | Formulation | Hydro | Beyond
---|---|---|---|---|---
| Source | catalog |  |  | GR
AMSS-NCKU [613, 614, 615, 552] | Y | – | BSSN/Z4c | – | Y
BAM [616, 582, 617, 618] | – | [619, 620, 524] | BSSN/Z4c | Y | –
BAMPS [621, 622, 623] | – | – | GHG | Y | –
COFFEE[624, 625] | Y | – | GCFE | – | Y
Dendro-GR [626, 627] | Y | – | BSSN/CCZ4 | – | Y
Einstein Toolkit [628, 629] | Y | – | BSSN/Z4c | Y | Y
∗Canuda [366, 367, 630] | Y | – | BSSN | – | Y
∗IllinoisGRMHD [631] | Y | – | BSSN | Y | –
∗LazEv [515, 632] | – | [633, 634, 635, 636] | BSSN/CCZ4 | – | –
∗Lean [637, 638] | Partially | – | BSSN | – | Y
∗MAYA [639] | – | [639] | BSSN | – | Y
∗NRPy+ [640] | Y | – | BSSN | Y | –
∗SphericalNR [641, 642] | – | – | spherical BSSN | Y | –
∗Spritz [643, 644] | Y | – | BSSN | Y | –
∗THC [645, 646, 647] | Y | [619] | BSSN/Z4c | Y | –
∗WhiskyMHD [648] | – | [649] | BSSN | Y | –
ExaHyPE [650] | Y | – | CCZ4 | Y | –
FIL[651] | – | – | BSSN/Z4c/CCZ4 | Y | –
GR-Athena++ [652] | Y | – | Z4c | Y | –
GRChombo [653, 654, 655] | Y | – | BSSN/CCZ4 | – | Y
HAD [656, 657, 658] | – | – | CCZ4 | Y | Y
Illinois GRMHD [659, 660] | – | – | BSSN | Y | –
MANGA/NRPy+ [661] | Partially | – | BSSN | Y | –
BH@H/NRPy+ [640, 662] | – | – | BSSN | – | –
MHDuet [663, 664, 665] | Y | – | CCZ4 | Y | Y
SACRA [666, 667, 668, 669, 670] | – | [671] | BSSN/Z4c | Y | Y
SACRA-SFS2D [672, 673] | – | – | BSSN/Z4c | Y | –
SpEC [674, 523] | – | [521, 523, 675] | GHG | Y | Y
SpECTRE [676, 677] | Y | – | GHG | Y | –
SPHINCS_BSSN [678] |  | – | BSSN | SPH | –
Table 7:  List of numerical relativity codes. We indicate if a code is open-
source, if it has been used to produce public gravitational waveform catalogs,
the formulation of Einstein’s equation used (GHG: generalized harmonic, BSSN:
Baumgarte-Shapiro-Shibata-Nakamura, CCZ4 / Z4c variants of the Z4 formulation,
GCFE: generalised conformal field equations), if a code implements general
relativistic hydrodynamics, and if it is capable of simulating compact
binaries beyond general relativity. An asterisk indicates codes that are
either (partially) based on the open-source Einstein Toolkit or are co-funded
by its grant. Note this table was created jointly for this paper and Ref.
[520].

In particular, LISA sources present a formidable challenge to NR (see Sec.
4.1.5), both with regard to the characteristic features of the target BBH
systems (e.g., high spins, large mass ratios, eccentric orbits, etc.) and the
quality of the GW data due to, e.g., high signal-to-noise ratios and sources
remaining in band for very long times. Computational simulations that meet the
accuracy requirements for these demanding BBH configurations, for sufficiently
long times, will require exascale computational resources. It will also
require continued research in NR to develop new algorithms, computational
techniques and software to tackle these challenges faced by existing code
bases. Some new codes are under development, and they are included in Table 7.

Given the complexity of NR codes, cross-validation is an important component
of code verification and it is by no means trivial that results obtained with
different code bases agree. They employ, for example, different theoretical
formulations of Einstein’s equations for the evolution and initial data and
different methods for the GW extraction. They also employ different numerical
techniques (e.g, pseudo-spectral or high-order finite differences
discretization) and numerical implementations, e.g., for the grid structure.
First systematic code comparisons and studies on the integration of NR with
data analysis and analytical approximations have been conducted in Refs. [679,
585, 592, 588, 587, 680]. For example, the NINJA project [585, 592] focused on
building a NR data analysis framework and to develop first injection studies.
The initial study presented a sample of $23$ NR waveforms of BBHs with
moderate mass ratios $q\lesssim 4$, considering only the dominant the
$\ell=|m|=2$ GW mode. The follow-up study increased to statistics to $60$ NR
waveforms. While no cross-validation between the NR waveforms were performed,
NINJA was crucial to identify technical and conceptual issues, that a
hybridization with PN or EOB may be needed and that more than only the
dominant GW mode may be necessary. The Samurai project [588] conducted a
detailed cross-validation of gravitational waveforms obtained with five
different NR codes. It focused on one BBH system, namely an equal-mass, non-
spinning, quasi-circular (eccentricity $\lesssim 0.0016$) binary completing
about six orbits before their merger. Focusing on the dominant $\ell=|m|=2$
multipole, it was found that the waveforms’ amplitude and phase agree within
numerical error. It was also found that these NR waveforms were
indistinguishable for SNR $\leq 14$ and would yield mismatches of $\lesssim
10^{-3}$ for binaries with $M\sim 60M_{\odot}$ (using the anticipated detector
noise curves at the time). The Numerical-Relativity–Analytical-Relativity
(NRAR) collaboration [587] made important strides towards constructing more
accurate inspiral-merger-ringdown waveforms by combining analytical and NR
computations, and by defining new accuracy standards that would be needed for
GW parameter estimation. The collaboration pushed cross-validation of NR codes
to new frontiers by considering quasi-circular BBHs (with initial eccentricity
$\leq 2\cdot 10^{-3}$) with both unequal mass (up to $q\leq 3$) and moderately
spinning BHs (up to $\chi_{i}\sim\pm 0.6$) completing about $20$ GW cycles
before merger. The study also included a binary of non-spinning BHs with
$q=10$, the highest mass ratio for quasi-circular binaries at the time. For
the first time, the seven involved NR codes performed $22$ targeted
simulations to meet the required accuracy standards. They employed the same
analysis code to estimate the uncertainties due to the numerical resolution
and waveform extractions. The simulations exhibited a relative amplitude error
in the $\ell=|m|=2$ multipole of $\lesssim 1\%$ and a cumulative phase error
of $\lesssim 0.25$rad. Finally, Ref. [680] presented a comparison of targeted
simulations for the first GW event, GW150914, concluding that the waveforms
were sufficiently accurate to effectively analyse LIGO data for comparable
events. A comparison between different time evolution formulations of
Einstein’s equations, namely BSSN and Z4c, was performed in Refs. [681, 552].

The cross-validation studies between different NR codes have been crucial for
building confidence in their results. As summarized above, they have been
conducted in limited regions of parameter space (e.g., low mass ratios, some
moderate spins, low eccentricity). As the demand on NR waveforms increases,
such as covering a larger region of parameter space, including eccentric or
spin-precessing BBHs, including higher multipoles (typically up to $\ell=8$)
that are excited by asymmetric systems and further improving their accuracy,
extended validation studies become important. With the technical and
computational developments, the numerical data becomes more and more sensitive
to mismatches, e.g., in the initial data, evolution or wave extraction.
Therefore, a careful study of the initial data and initial parameters, of wave
extraction techniques (e.g., extrapolation of waveforms extracted at finite
radii vs. Cauchy characteristic extraction vs. Cauchy characteristic matching)
or the choice of asymptotic Bondi-van der Burg-Metzner-Sachs frame will be
needed [682].

The starting point for any simulation is the construction of constraint
satisfying initial data. After decades of work (see, e.g., Refs. [683, 529,
684, 531, 685, 686, 687, 533, 536] and references therein), multiple codes are
now capable of generating BBH initial data where the BHs are on a quasi-
circular orbit or where the orbit has non-zero eccentricity.

Several groups worldwide have created public catalogs of BBH merger
simulations with approximately $5,700$ waveforms available at the time of this
writing [521, 639, 633, 523, 634, 635, 636, 524]; see also Table 7. Nowadays,
NR waveforms survey several configurations in the extensive BBH parameter
space spanned by their mass ratios, spin magnitudes and directions, and
eccentricity. Fig. 3 illustrates the portion of the parameter space covered by
the publicly available waveforms at the time of this paper’s publication.
Focussing on quasi-circular binaries, it shows that the parameter space is
best sampled for moderate mass ratios and spins.

Figure 3:  Parameter space coverage of public waveform catalogs in the quasi-
circular limit. Shown are mass-ratio $q$, projections of the BH spins onto the
orbital angular momentum ($\chi_{1/2\,z}$), and magnitude of spin-components
orthogonal to the orbital angular momentum $(\chi_{1/2\,\perp}$) for waveform
catalogs of different NR groups. Data as of June 2023.

Consequently, there are three dimensions in the BBH parameter space — spin,
mass ratio and eccentricity — that require further waveform development. For
example, most runs with spinning BHs concentrate on spins $\chi\lesssim 0.9$,
because simulating BHs with spin closer to the extremal Kerr limit is
technically challenging. It requires a much higher numerical resolution and,
more importantly, the initial data is more difficult to construct. There are
some codes capable of constructing initial data for BHs with spins above
$\chi=0.95$; see, e.g., Refs. [533, 688, 689]. Another challenging regime in
BBH simulations is that of unequal mass binaries. Most simulations in the
public catalogs cover the regime $q\lesssim 8$ for spinning BHs, and up to
$q\lesssim 18$ for nonspinning binaries, as can be seen in Fig. 3. However,
simulations above $q=20$ are scarce: the first high-mass ratio simulations up
to $q=100$ [690] (see also [626]) and, more recently, up to $q=1024$ [691]
were achieved for head-on collisions. The frontier for quasi-circular
inspirals are mass ratios of $q=128$ [692, 693]. New techniques such as
worldtube excision to model intermediate mass ratios of $q=100$-$10^{4}$ are
under development [694, 695].

Most NR simulations to-date have focused on quasi-circular inspirals because
it is expected that any eccentricity present during a binary’s formation will
have been emitted through GW radiation when it close to merger [192]. However,
there are several astrophysical scenarios that predict eccentric binaries
[696, 697, 698, 699, 700, 701, 702, 703]. The NR community is beginning to
explore the eccentric parameter space [704, 705, 706, 707, 708, 709, 710, 711,
712, 713, 714, 715, 596, 716, 717, 636] as well as hyperbolic encounters [718,
719, 720, 721].

Simulating eccentric binaries is challenging for a few reasons: (i) unlike in
Newtonian gravity, there is no unique definition of eccentricity, and NR
groups have explored different definitions [722, 711, 723, 724, 725]; (ii)
increasing the eccentricity to for a fixed separation requires that the BBH
would have to start further and further apart to undergo several GW cycles
before merger which yields higher numerical cost. Despite these challenges,
eccentric NR runs are becoming increasingly available, although the parameter
space in mass ratio and spin is not (yet) well sampled, particularly for
precessing systems.

In summary, the BBH parameter space is currently covered as follows:

  * $\bullet$

Non-spinning BBH mergers with mass ratio $q\leq 18$ [639, 633, 634, 521, 523,
586, 726, 587] — see [727, 692] for binaries with mass ratios $q=64,100$, and
$128$.

  * $\bullet$

BBHs with moderate random spins and $q\leq 8$ [639, 521, 633, 634, 587, 524],
aligned-spin binaries of $\chi\leq 0.85$ with $q\leq 18$ [522, 728], and a few
with very high spins ($\chi\sim 0.99$) [689, 583, 729].

  * $\bullet$

Eccentric BBHs with $q\leq 10$ [704, 708, 710, 596].

Note that BBH simulations in GR are scale-invariant with respect to the total
mass of the system; therefore, a BBH waveform can be re-scaled for both
ground-based and space-based gravitational wave detectors without having to
repeat the NR computation. However, the mass scaling will determine the
physical starting frequency of the signal, and therefore whether the NR
waveform covers the entire sensitivity band of the detector.

####  4.1.4 Environmental Effects

The NR community has made great strides in simulating binaries in a matter-
rich environment. This Whitepaper focuses on ways the environment directly
impacts the predicted gravitational waves; for EM signatures of BBH mergers,
we refer to the Astrophysics Working Group White Paper [3] and references
therein.

The gravitational waves emitted during the merger of a BBH in an accretion
disk are likely unaffected by the accreting matter [730], although Ref. [731]
indicates that accretion disk densities greater than $10^{6}-10^{7}$ g cm-3
would alter the coalescence dynamics enough to be relevant for ground-based GW
detectors.

In addition to EM counterparts, there are other possible sources for LISA
modeled by NR. One such class of signals arises from the magnetorotational
collapse of a supermassive stars; see, e.g., [732] for a recent review and
[733] for early NR work. Recently the GW signatures from such collapsing
supermassive stars were calculated in [734, 735], where it was shown that LISA
could observe such events out to redshift $z\simeq 3$. Ref [735] also computed
possible EM counterparts and predicted that these systems could be sources of
very long gamma-ray bursts.

Another potential source of gravitational waves detectable by LISA that is
currently unmodeled comes from instabilities in accretion disks around BHs.
The particular case of the Papaloizou-Pringle instability [736] of self-
gravitating disks has been studied in [737] without BH spin, and in [738]
including spin. In [738] it was shown that BH spin could potentially increase
the duration of the near monochromatic signal from the instability, and that
LISA could detect gravitational waves generated by the Papaloizou-Pringle
instability around $10^{5}M_{\odot}$ BHs out to redshift $z\simeq 1$.

Finally, for a detailed discussion of environmental effects and the challenges
they could present for precision tests of GR, such as dark matter
environments, we refer the interested reader to the Fundamental Physics
Whitepaper [323] and references therein.

####  4.1.5 Challenges in NR

Improvements in hardware and numerical techniques continue to speed up
simulations, yet some corners of the binary parameter space continue to demand
challenging simulations that are costly in terms of runtime and computational
resources. Three of the most pressing challenges facing NR are i) producing
waveforms that cover the anticipated parameter space of unequal mass, highly
spinning, and eccentric binaries, ii) producing long-lived waveforms with
sufficient numbers of GW cycles, and iii) producing these waveforms at
standards of accuracy set by LISA’s sensitivity.

##### Parameter coverage:

The BBH parameter space for the LISA mission, as summarized in Table (2) has
mass ratios that span from 1:1 to 1:1000, BH spins from 0 to 0.998, and a wide
range of eccentricities. As discussed in $\S$ 4.1.3, almost all NR simulations
to date are for $q\leq 8$, $\chi\leq 0.8$ and $e\approx 0$; and, therefore,
there is significant work to be done to supply waveforms that cover the full
potential parameter space for LISA. We also discuss LISA waveform catalogs.

Producing waveforms for systems with large mass ratios is a challenge for the
broader waveform community with various approaches existing to bridge the gap
between NR solutions and small mass ratio approximations. This is especially
relevant for IMRIs. High-mass ratios are demanding due to the need to resolve
the smaller mass black hole and provide appropriate gauge conditions [693]. As
discussed in $\S$ 4.1.3, several new codes are being developed by the NR
community with the goal of having simulations with large mass ratios ($q>50$)
be routinely possible. In addition to pushing the NR capacity to larger mass
ratios, new methods for modeling IMRIs are currently being explored [739,
694], which combine black-hole perturbation theory and NR techniques to
significantly increase the numerical efficiency of simulations. Coordination
between the GSF and NR communities as well as the construction of open source
platforms to share NR [523, 633, 639] and GSF waveforms, such as the Einstein
Toolkit [628] and the Black Hole Perturbation Toolkit [740], will streamline
and accelerate such endeavors.

A second challenge is producing many generic simulations with BH spins that
are very large in magnitude, i.e., approaching 1. The simplest approach to
solving the initial data problem is to use the Bowen and York method [741,
529]; however, this method cannot construct BHs with high spins, $\chi\gtrsim
0.93$ [742, 743, 744]. In addition, resolving the region of spacetime near the
horizons requires computationally expensive, very high resolution. This is
exacerbated in evolution methods using excision to handle the singularities
[583]. As mentioned previously, NR codes have successfully achieved some high
spins, including several aligned-spin binaries with spins ($\chi\sim 0.99$)
[689, 583, 729] by using new formulations that move beyond the conformally
flat ansatz for puncture methods [536] and new techniques for handling the
excision region [583].

Finally, several astrophysical scenarios produce non-zero eccentricity in the
LISA frequency band [745]. NR is capable of and has been producing eccentric
runs [704, 708, 710, 596, 717, 746]; the challenge is simulations must start
with the BHs more widely separated than in quasi-circular configurations. This
is again a computational cost.

The source modeling community has continued to improve the modeling of these
sources, developing more accurate inspiral-only models [747, 748, 749, 750,
751], and inspiral-merger-ringdown models that combine analytical
approximations [752] with eccentric NR waveforms [704, 709, 710, 711, 712,
713, 714, 715, 596].

##### Gravitational wave cycles:

Even after new NR codes attain the ability to simulate accurate, generic high
mass ratios, advances in approximate methods will still be needed. LISA will
observe BBH mergers for months or even years, and it is unrealistic (and
unnecessary given the accuracy of approximate waveform models calibrated with
NR [753, 754, 728, 589, 522]) to use NR to describe the entire evolution of
binary systems, from inspiral to ringdown. Semi-analytic models of the
inspiral-merger-ringdown have been very effective in ground-based GW
astronomy, notably Effective-One-Body [755, 756, 757, 758, 759, 760],
Phenomenological [753, 728, 761, 762, 763] and surrogate models [764, 590,
765, 766]. Additionally, combining GSF methods to model the inspiral evolution
with NR for the late-inspiral, merger and ringdown stages as described in
[694] is a promising avenue of study.

##### Accuracy

Massive BBH mergers in LISA could have high (on the order of 1000) signal-to-
noise (SNR) ratios and NR must model these loud signals accurately enough not
only to infer the correct source properties, but also to subtract the signal
from LISA data, revealing quieter, high-interest signals that may lay
underneath. The full implications of these requirements for NR waveform
accuracy have yet to be assessed, but an initial study [2] suggests that NR
waveforms will require a substantial increase in accuracy, compared to today’s
most accurate waveforms, to be indistinguishable from a high-SNR observation
with the same source properties. A similar assessment for next-generation
ground-based detectors concluded that NR waveforms will need an order of
magnitude more accuracy to avoid bias in the inferred source properties for
the high signal-to-noise sources that these detectors will observe [1]. The
impact of using insufficiently accurate templates is highlighted in Figure 4.

Figure 4:  Illustration of how numerical errors in a template may leave a
residual in high-SNR signals. The left panel represents a BBH at high SNR
(detector mass $M_{\mathrm{d}et}=5\times 10^{6}M_{\odot}$, aligned spin of
$\chi=0.2$ on the larger BH at distance $D_{L}=30\,$Gpc and at inclination
$\iota=15^{\circ}$). The black line is the characteristic strain
($2f\tilde{h}(f)$) obtained at high numerical resolution from all spherical
harmonic modes. The blue dashed line is the corresponding quantity obtained
from a low resolution simulation, with the light blue line denoting the
difference to the high-resolution waveform. The solid red line is a high
resolution template containing only the $(l,m)=(2,2)$ mode, with the light red
line showing the difference to the black line. If this system were searched
for with the low-resolution (or (2,2)-mode only) waveforms, then the shaded
lines would remain as residuals after subtraction. The SNR of the residual due
to resolution in the left panel is 342, clearly indicating the need for higher
accuracy NR simulations. The right panel represents a more frequently expected
system, where the SNR is 200 and $\iota=0^{\circ}$. Here, the residual due to
resolution has an SNR of 2.4. Figure adapted from [2].

These accuracy requirements are even more challenging for NR if the BHs are
precessing, if one or both BHs are spinning nearly as rapidly as possible, or
if one BH is much more massive than the other. Let us illustrate how quickly
computational cost rises, by quantifying computational cost increases relative
to a fiducial baseline simulation with mass-ratio $q_{0}$ and initial
frequency $M\Omega_{0}$. Cost increases arise for multiple reasons.

First, the explicit evolution schemes used in NR codes are subject to the
Courant-Friedrichs-Levy (CFL) condition, which limits the possible time step
to the dynamic time of the smaller BH, $\Delta t\propto m_{2}$. As the mass
ratio increases, the number of time-steps per unit $M=m_{1}+m_{2}$ of
evolution time will therefore increase as

$f_{\mathrm{C}FL}\sim\frac{q}{q_{0}}.$ (7)

Second, the total inspiral time of a BBH can be estimated by PN expansions
[767] as $T\propto 1/\nu(M\Omega)^{-8/3}$, where $M\Omega$ is the initial
dimensionless frequency of a binary. A change in mass ratio and initial
frequency will modify the computational cost in proportion to $T$, i.e. by a
factor

$f_{\Omega_{\mathrm{l}ow}}=\left(\frac{\nu}{\nu_{0}}\right)^{-1}\left(\frac{\Omega}{\Omega_{0}}\right)^{-8/3}.$
(8)

Thirdly, LISA has higher accuracy requirements than current NR simulations. We
denote by $f_{\tau}(A)$ the increase in computational cost that is needed to
reduce the truncation error $\tau$ by a factor $A>1$, i.e.
$\tau\rightarrow\tau/A$.

A fourth cause for increased computational cost arises from the desire of
achieving a certain overall phase error, which is important for GW waveform
modeling. For longer simulations, phase accuracy must be preserved over more
cycles. From PN theory, the orbital phase to merger scales as
$\Phi\propto\nu^{-1}(M\Omega_{i})^{-5/3}$. To achieve a fixed overall phase-
error (say, 1 radian), the relative phase error must decrease $\propto
1/\Phi$, yielding an increase in computational cost by a factor

$f_{\Phi}\sim
f_{\tau}\left(\left(\frac{\nu}{\nu_{0}}\right)^{-1}\left(\frac{\Omega}{\Omega_{0}}\right)^{-5/3}\right).$
(9)

The challenge now is that in order to achieve higher mass ratio simulations
starting at similar (or lower) initial frequency, and possibly also achieving
higher accuracy, each of the factors just outlined increases very rapidly and
that these factors multiply for an overall cost increase by a factor

$C=f_{\mathrm{C}FL}\,f_{\Omega_{\mathrm{l}ow}}\,f_{\tau}\,f_{\Phi}.$ (10)

To obtain concrete numerical estimates, we assume that $f_{\tau}$ is a power
law, $f_{\tau}(A)=A^{\alpha}$. This assumption is satisfied for finite
difference codes, where a convergence order of $k$ in $(3+1)$ dimension yields
$\alpha=4/k$. In terms of $\alpha$, the overall cost increase is given by

$C\sim\left(\frac{q}{q_{0}}\right)^{2+\alpha}\,\left(\frac{\Omega_{0}}{\Omega}\right)^{8/3+5\alpha/3}\,A^{\alpha},$
(11)

where we have also used $\nu\approx 1/q$, which is valid at large $q$.

We can now make concrete estimates: Increasing the mass-ratio by a factor of
$2$, lowering the initial frequency by a factor of $1/2$, and increasing
accuracy by a factor of $2$, and assuming the NR code under consideration
maintains perfect 8-th order convergence into the regime of longer, more
accurate simulations at higher mass-ratio (i.e. $\alpha=1/2$), increases the
computational cost two orders of magnitude. More ambitiously, increasing the
mass-ratio by a factor of 10, reducing initial frequency by a factor $1/4$ and
increasing accuracy by a factor $10$ increases the computational cost by five
orders of magnitude. Codes based on spectral methods, e.g. SpEC, will likely
result in a smaller $f_{\tau}$; even so, increases in computational cost are
very challenging as one goes to longer simulations at higher mass-ratio.

The NR codes that have been instrumental in providing public BBH waveforms use
either a spectral method or box-in-box mesh refinement. These methods scale
well to the order of 1000 cores, but they cannot scale to hundreds of
thousands or millions of cores. For example, codes using high-order finite-
difference stencils incur ever-increasing inter-core communication costs as
the number of cores increases and this cost eventually becomes prohibitive.
These codes are parallelization bound and do not parallelize beyond current
usage. Since individual computer cores are no longer dramatically increasing
in speed (Moore’s law has moved from single core performance into increasing
the number of cores), it therefore seems unlikely that these codes will be
able to simulate long BBH inspirals at mass-ratios significantly larger than
$10$. In the section on next-generation codes, we explore how this is being
improved.

For a more general discussion on waveform accuracy, see $\S$ 3.2.

##### NR Catalogs

The LISA data analysis groups could have access to NR waveforms via a bespoke
LISA catalog or as a public catalog. A bespoke catalog would have NR waveforms
created specifically to meet those aforementioned requirements, potentially
interfacing and merging existing catalogs and extending with its own dedicated
entries (similar to the LVK catalog [591]). If a bespoke LISA catalog is not
created or deemed unnecessary due to a possible future abundance of NR public
waveforms (current catalogs are indicated in Table 7) we will need a software
interface with existing NR waveform catalogs that is curated and cultivated
with time. A full catalog of NR waveforms achieving all requirements will take
time, possibly until the launch of LISA. We put forth the following priorities
for the NR waveforms: 1) an accuracy assessment of current and near-term
achievable NR waveforms, 2) the assessment of current and near-term accessible
parameter space, 3) the development of NR waveforms for the most likely LISA
events at the required accuracy, 4) assessment of systematic bias in simulated
LISA signals with currently available waveforms, and 5) the development of NR
waveform for all potential parameters at the necessary accuracy.

###  4.2 Weak field approximations (post-Newtonian/post-Minkowskian)

Coordinator: Laura Bernard and Chris Kavanagh
Contributors: A. Antonelli, G. Faye, J. Garcia-Bellido, M. Haney, F.
Larrouturou, A. Le Tiec, O. Long

####  4.2.1 Description

The PN formalism is an approximation method in GR that is well suited to
describe the orbital motion of—and the GW emission from—binary systems of
compact objects, in a regime where the typical orbital velocities are small
compared to the vacuum speed of light $c$ and the gravitational fields are
weak. It has played and continues to play a central role in the construction
of template banks for the detection and analysis of GW signals generated by
binary systems of black hole and/or neutron star, which are now routinely
observed by the LVK collaboration’s detectors [193, 768, 769, 770, 771, 210,
114]. More precisely, PN results are at the core of several classes of
template banks, including phenomenological waveforms and effective one-body
waveforms (see Secs. 4.5 and 4.6).

In PN theory, the general relativistic corrections to the Newtonian motion and
to the leading GW emission (i.e. Einstein’s quadrupole formulas) are computed
in a systematic manner in powers of the small PN parameter $v^{2}/c^{2}\sim
GM/(c^{2}r)$, where $v$ and $r$ are the typical binary relative velocity and
separation, $M$ is the sum of the component masses, and $v^{2}\sim GM/r$ for
bound orbits. Indeed, the most promising sources of gravitational waves for
existing and planned interferometric detectors are bound systems of compact
objects. By convention, a contribution of “$n$PN” order refers to equation-of-
motion terms that are $O(1/c^{2n})$ smaller than the Newtonian acceleration,
or, in the radiation field, smaller by that factor relative to the standard
quadrupolar field.

The PN approximation dates back to the pioneering works of Lorentz and Droste
[772], as well as Einstein, Infeld and Hoffman [773], who computed the leading
1PN corrections to the Newtonian equations of motion in a system of $N$ point
masses. During the 1980s, those results were extended to 2.5PN to provide a
rigorous basis for interpreting binary pulsar observations [774, 775]. For
binary systems of compact objects, the state of the art corresponds to
gravitational waveforms that include all of the relativistic corrections up to
4PN order, in the simplest case of nonspinning bodies moving along a sequence
of quasi-circular orbits.

Another complementary weak field approximation can be found by relaxing the
small velocity assumption of the PN approximation and demanding only that
$GM/(c^{2}r)\ll 1$. This scheme, known as the PM approximation, has a long and
venerable history, in the context of both unbound scattering where velocities
can be arbitrarily large [776, 777] and far-zone–near-zone matching for bound
orbits within certain PN schemes (see e.g., Ref. [767]).

As discussed below, much work has recently been dedicated to (i) push this
accuracy to 4.5PN order, (ii) include the effects of the spins of the compact
objects (see also Table 8), (iii) generalize the results from circular to
eccentric orbits, and (iv) explore the overlap and synergies with other
approximation methods (PM approximation, small-mass-ratio approximation,
effective one-body model). The reader is referred to the review articles [778,
779, 780, 767, 781, 782, 783, 784, 785] and to the texbooks [786, 787] for
more information.

| Dynamics | Dissipative flux
---|---|---
PN order | non-spinning | spinning | non-spinning | spinning
|  | SO | SS | higher spins |  | SO | SS | higher spins
0 | $\surd$ | - | - | - | - | - | - | -
1 | $\surd$ | - | - | - | - | - | - | -
1.5 | - | $\surd$ | - | - | - | - | - | -
2 | $\surd$ | - | $\surd$ | - | - | - | - | -
2.5 | $\surd$ | $\surd$ | - | - | $\surd$ | - | - | -
3 | $\surd$ | - | $\surd$ | - | - | - | - | -
3.5 | $\surd$ | $\surd$ | - | $\surd\ (S^{3})$ | $\surd$ | - | - | -
4 | $\surd$ | - | $\surd$ | $\surd\ (S^{4})$ | $\surd$ | $\surd$ | - | -
4.5 | $\ast$ | $\surd$ | - | $\surd\ (S^{3})$ | $\surd$ | - | $\surd$ | -
5 | $\ast$ | - | $\surd$ | $\surd\ (S^{4})$ | $\surd$ | $\surd$ | - | -
5.5 | $\ast$ |  |  | $\surd\ (S^{5})$ | $\surd$ | $\surd$ | $\surd$ | -
6 |  |  |  | $\surd\ (S^{6})$ | $\surd$ | $\surd$ | $\surd$ | $\surd\ (S^{3})$
6.5 |  |  |  | $\star$ | $\surd$ | $\surd$ | $\surd$ |
7 |  |  |  | $\star$ | $\surd$ |  |  |
Table 8: State-of-the-art of known PN results for both the conservative and
dissipative dynamics as well as for the gravitational flux. Contrary to the
text, everything is stated as absolute order. For example, the 6PN absolute
order for the non-spinning flux in the table corresponds to the 3.5PN relative
order results as stated in the text. Hereditary effects are only known as non-
local contributions. An instantaneous expression can be obtained by performing
a low-eccentricity expansion. SO and SS refer to spin-orbit and spin-spin
interactions respectively.
$\ast$ means that only a partial result is known at those orders. At 5 and
5.5PN order, the results were obtained by the combination of PN traditional
techniques with scattering amplitudes and self-force.
$\star$ means that the dynamics is known at all leading order in the spin.

####  4.2.2 Suitable for what sources?

The PN approximation is well suited to describe the dynamics and gravitational
emission of compact binary systems in the early inspiral stage. While it can
in principle be applied to model BH binaries of any mass ratio, it is expected
to be more accurate for comparable to intermediate mass BHBs, as the
gravitational field should remain small until the final cycles. For mass ratio
above $q\gtrsim 100$, when the small body remains for long timescales close to
the supermassive BH, the gravitational field is strong and the weak-field
approximation will lose accuracy (see Fig. 2).

The PM approximation is best suited to close hyperbolic encounters (CHE),
which could be detected with LISA for a population of massive PBH in dense
clusters forming part of the halos of galaxies, as studied in Refs. [788,
789]. It also provides valuable resummations of terms in the PN approximation
and is used to improve bound compact binary coalescence models when combined
with PN results.

PN waveforms are ideally suited to model binaries that are in an early stage
of their inspiral while in the LISA band such as SOBHBs and GBs.

####  4.2.3 Status

#####  4.2.3.1 Equations of motion and waveforms without spin.

For an isolated, non-spinning compact binary system, the phase evolution and
GW modes during the inspiral stage are currently respectively known up to
4.5PN and 4PN [790, 791, 792, 793, 794, 795, 796, 797, 798]. As for the LISA
requirements, it was established in [799] that the 2PN waveform will be enough
for the parameter estimation of roughly 90% of the detectable SOBHBs, while
the 3PN waveform phasing is needed for the systems that will merge within the
mission lifetime. Nevertheless it will be quite valuable to push the current
accuracy up to 4.5PN order for at least two main reasons. The first one is
that it may be necessary for the parameter estimation of black hole binaries
with large masses. The second one is that this accuracy will be very
beneficial for the calibration of numerical waveforms, and for testing the
second-order self-force computations.

For quasi-circular orbits, the gravitational phase $\phi(t)$ is computed via
the so-called energy _flux-balance equation_ ,
$\frac{\mathrm{d}E}{\mathrm{d}t}=-\mathcal{F}_{\text{GW}}$, where $E$ is the
conserved energy and $\mathcal{F}_{\text{GW}}$ is the flux of radiation
emitted. Both quantities are functions of the gravitational phase, which is
thus obtained by a simple integration. The amplitude and polarizations are
obtained by the same radiative multipole moments that are used when computing
the flux, as explained later. For recent reviews, see [800, 767, 801, 781,
783].

At 4PN precision, three different techniques are used to compute the binding
energy: the canonical Hamiltonian formulation of General Relativity [802, 803,
804, 805, 806], the Fokker Lagrangian approach [807, 808, 809, 810, 811] and
effective field theory (EFT) [812, 813, 814, 783, 815, 816, 817]. All of these
methods derive the energy as a conserved quantity associated with the
equations of motion for two point-like particles, and naturally give
physically equivalent results. The Hamiltonian approach has been pushed up to
6PN, yielding incomplete results [818], see Sec. 4.2.3.4. Exploiting synergies
between traditional PN methods and EFT, the logarithmic contributions in the
energy have also been computed up to 7PN [819].

The flux of gravitational radiation can be expressed as a generalisation of
the famous Einstein quadrupole formula [820, 821], as [822]

$\displaystyle\mathcal{F}_{\text{GW}}$ $\displaystyle=\sum_{\ell\geq
2}\frac{G}{c^{2\ell+1}}\left[\alpha_{\ell}\,U_{L}^{(1)}U_{L}^{(1)}+\frac{\beta_{\ell}}{c^{2}}\,V_{L}^{(1)}V_{L}^{(1)}\right]$
$\displaystyle=\sum_{\ell\geq
2}\frac{G}{c^{2\ell+1}}\left[\alpha_{\ell}\,I_{L}^{(\ell+1)}I_{L}^{(\ell+1)}+\frac{\beta_{\ell}}{c^{2}}\,J_{L}^{(\ell+1)}J_{L}^{(\ell+1)}\right]+\mathcal{F}_{\text{tails}}\,,$
(12)

where the $U_{L}$ (resp. $V_{L}$) are the mass (current) radiative multipole
moments and the $I_{L}$ (resp. $J_{L}$) are the mass (current) source
multipole moments, the $\\{\alpha_{\ell},\beta_{\ell}\\}$ are collections of
numbers and the superscript $(n)$ denotes $n$ time derivatives. To connect
radiative moments to sources moments, their non-linear interactions (occurring
during the propagation from the source to the detector) have been singled out
in $\mathcal{F}_{\text{tails}}$. The radiative moments are also directly
involved in the computation of polarisations [822], and thus their
determination is crucial.

The contribution $\mathcal{F}_{\text{tails}}$ is known up to 4.5PN order,
using traditional PN methods [823], and was confirmed by an independent PN re-
expansion of resumed waveforms [824]. The computation of the required moments
is currently done using the Multipolar-post-Minkowskian-post-Newtonian (MPM-
PN) algorithm [767]. The major piece of this work is the derivation of the
mass quadrupole ($\ell=2$) up to 4PN order: the main result has been obtained
[790, 792, 793, 794, 795], including the tail-of-memory contribution [825].
The other moments involved are the mass octupole and the current quadrupole,
already known at the required order (3PN) [826, 791], and higher moments that
are either known or trivial (i.e. needed at the Newtonian order only), see
notably [827]. Collecting all these results, the complete flux of
gravitational radiation at infinity is fully known at 4PN [798] (and 4.5PN in
the case of circular orbits [823]). In the case of horizonless objects such as
neutron stars, the flux-balance equation is used to compute the gravitational
phase for circular orbits at 4.5PN [797]. However, if at least one object is a
black hole, one must include an extra contribution due to the absorption of
gravitational waves by the black hole horizon(s) [828, 829, 830, 831, 832].
Note that two other methods are able to deal with high PN multipole moments,
both currently developed up to 2PN: the direct integration of the relaxed
equations (DIRE) [833] and EFT [834]. As half-PN orders are often easier to
compute than integer ones, one can probably push the results up to 4.5PN
order, including the amplitudes of the gravitational modes. A first step in
this program is to control the dissipative, radiation reaction effects in the
equations of motion. Table 8 summarizes the current knowledge regarding PN
dynamics and flux without spin.

Finally, there have been some works trying to map the BMS (Bondi-Metzner-
Sachs) flux-balance laws to equivalent PN results in the harmonic gauge [835].
By transforming the metric from harmonic coordinates to radiative Newman-Unti
(NU) coordinates, they were able to obtain the mass and angular momentum
aspects and the Bondi shear as a function of the quadrupole moment [836]. In
particular, they rederive the displacement memory effect (see also [564]) and
provide expressions for all Bondi aspects relevant to the study of leading and
subleading memory effects [837].

#####  4.2.3.2 Equations of motion and waveforms with spin.

In the last decade, the effect of the proper rotation of the components of
binary systems have been included in the orbital equations of motion to 4PN
order, extending previous works dating back to the 70’s [838]. In the post-
Newtonian terminology, the physical spin is rescaled as $S/c$. The leading
spin interactions, the so-called spin-orbit interaction (SO), couple the mass
$m_{1}$ of a particle with the spin of the other, hence those contributions
scale as
$\frac{Gm_{1}}{c^{2}}\times\frac{S_{2}}{c}\sim\mathcal{O}(\frac{1}{c^{3}})$
and are regarded as being of 1.5PN order. Meanwhile the first effect due to
spin-spin interaction, scaling as
$\frac{G}{c^{2}}\times\frac{S_{1}}{c}\times\frac{S_{2}}{c}\sim\mathcal{O}(\frac{1}{c^{4}})$,
is quadratic in spin (SS) and arises at 2PN order. We refer to the successive
subleading PN contributions to a given spin interaction (e.g. SO or SS) as the
next-to-leading (NL), next-to-next-to-leading (NNL) and so on. Moving beyond
the 2PN two-body dynamics including SO interactions [839], it becomes
increasingly difficult to resort to an explicit model, where the bodies are
described as small balls of fluid [840]. The most efficient strategy consists
in adopting an effective point of view and considering that the two objects
are point particles endowed with a classical spin.

The definition of spin in GR was first introduced by Mathisson [841], and
rephrased later by Tulczyjew [842]. The covariant equations of evolution
obeyed by a spinning test particle were obtained in their usual form by
Papapetrou [843] and later generalized by Dixon to bodies endowed with higher-
order multipole moments [844]. The equations of motion and precession served
as a starting point to investigate the dynamics of binary systems at next-to-
leading (NL: 2.5PN) [845, 846] and NNL (3.5PN) order linear in spin (SO
interactions) [847], assuming that all Dixon moments other than the masses and
the spins vanish. The metric is obtained by solving iteratively Einstein’s
equations in harmonic coordinates [848] for the pole-dipole stress-energy
tensor. Later, to compute the SS corrections at NL (3PN) order [849], the
contribution of the Dixon quadrupoles was added, so as to account for the self
deformation of the bodies produced by their own spins.

The direct computation of a generalized Lagrangian (i.e. depending on the
accelerations) for two spinning particles, using EFT in harmonic coordinates,
was performed in stages from the mid-2000’s on by integrating out the
gravitational field entering the full ‘field plus matter’ Lagrangian with the
help of standard Feymann diagram expansions [850, 851, 852]. The NL
$S_{1}$–$S_{2}$ (3PN) interactions [853, 854] were computed first, before the
NL $S_{1}$–$S_{1}$ (3PN) contributions [855] and NL SO (2.5PN) interactions
[856, 857] were considered. Equivalent results were obtained in parallel by
means of Hamiltonian methods based of the Arnowitt-Deser-Misner (ADM)
formulation of GR [858, 859, 860].

After the canonical treatment of the spin was better understood [860, 861],
the computations of the Hamiltonian were pushed up to NNL (3.5PN) order for
the SO [862, 863], to NNL (4PN) order for the $S_{1}$–$S_{2}$ interactions
[858, 864] and to NL (3PN) order for the $S^{2}$ interaction [865]. A fully
equivalent Lagrangian for the latter effects was found using EFT [866]. Once
the EFT degrees of freedom and the gauge choice corresponding to the spin-
supplementary condition were clearly identified, the PN leading terms to all
orders in spins were computed [867, 868, 869, 870] and the spin part of the
harmonic Lagrangian was completed up to the NNNL (5PN) order [871, 872].

Note that reaction forces, dissipative in essence, are absent from the
previous treatments, which exclusively describe the conservative dynamics.
Nonetheless, they can be computed from balance equations [873], or more
directly from the Papapetrou evolution equations, or from the conservation of
an appropriate stress-energy tensor [874]. One may even construct a Schwinger-
Keldysh Lagrangian where each degree of freedom is formally doubled. This was
done for the NL SO and SS effects at 4PN and 4.5PN respectively in Refs. [875,
876], and to all orders in spin in Ref. [877] (see also the approach developed
for ADM Hamiltonians [878] at leading SO and SS orders).

Knowing the near-zone dynamics of two spinning particles [848], one can insert
the corresponding PN metric into the right-hand side of Einstein’s equations
in harmonic coordinates, which yields the expression of the effective non-
linear source entering the integrand of the multipole moments, as defined in
the multipolar post-Minkowskian (MPM) formalism [879], or, equivalently, in
EFT [880]; hence one gets the relevant mass, current and gauge moments in the
spinning case [881]. The expression for the radiative moments in terms of the
former quantities has been derived for general isolated systems [826], but the
hereditary terms therein are more delicate to evaluate than for mere point-
mass particles on circular orbits, since the orbital plane is now generically
precessing. They have however been handled for the leading [882] and NL [883]
SO interactions, which required the full integration of the (approximately)
conservative dynamics up to the considered order.

Finally concerning the gravitational radiation, the GW flux has been obtained
to NNL (4PN) order both for the SO and SS contributions [849, 884]. The GW
phase, built from the radiative moments and the Noetherian energy, has been
obtained within the MPM formalism at the NNL (4PN) order for the SO
contributions [885, 883] (in continuation of previous works at leading order
(2PN) [886, 887] and NL (3PN) orders [888]), at the NNL (4PN) order for the SS
terms [849], and at leading order for the spin cube contributions [889].
Similar results were obtained with the EFT formalism [890] for which state-of-
the-art results are also the NNL (4PN) order [884]. The waveform modes have
been obtained to NNL (3.5PN) order [891] in the quasi-circular orbit
approximations. By contrast, the radiation amplitude has been computed up to
the NL (2PN) order only [892], although it is currently provided in a ready-
to-use form for precessing quasi-circular orbits to the even lower 1.5PN order
[893].

See Table 8 for a summary of on the current knowledge about the PN spinning
dynamics and flux.

#####  4.2.3.3 Eccentric-orbit waveforms.

The modeling of inspiral waveforms from compact binaries in eccentric orbits
commonly relies on quasi-Keplerian parametrization (QKP) as a semi-analytic
representation of the perturbative, slowly precessing post-Newtonian motion.
The conservative motion as well as instantaneous and hereditary contributions
to the secular (“orbital averaged”) evolution of the orbital elements are
known to 3PN order, in both modified harmonic and ADM-type coordinates [894,
895, 896, 897, 898] and to 4PN order in the ADM-type coordinates only [899].

The complete 3PN-accurate GW amplitudes from non-spinning eccentric binaries
have been derived using the multipolar post-Minkowskian formalism, including
all instantaneous, tail and non-linear contributions to the spherical harmonic
modes [900, 901, 751]. Going beyond the usual approximation of radiation
reaction as an adiabatic process (and the associated “orbital averaging” in
QKP), [894, 902, 903, 901] provide post-adiabatic, oscillatory corrections to
the secular evolution of GW phase and amplitude.

For spinning eccentric binaries that have component spins aligned with the
orbital angular momentum, the effects of spin-orbit and spin-spin couplings on
the binary evolution and gravitational radiation have been worked out to
leading order in QKP (i.e., up to 1.5PN and 2PN in the equations of motion)
[904]. Subsequent work has aimed to extend the treatment of spins in QKP to
higher PN orders [905, 906]. The waveforms modes to NL (3PN) order for aligned
spins and eccentric orbits have been derived in [757, 907], including tail and
memory contributions. While the instantaneous contributions were derived for
generic motion, the hereditary contributions were computed first in a small-
eccentricity expansion, then extended to larger eccentricities using a
resummation. Regarding precessing eccentric systems, gravitational waveforms
have been obtained to leading order in the precessing equation [908]. A fully
analytical treatment has also been proposed up to 2PN order in the spin,
including higher harmonics [909].

In practice, the semi-analytic approach of QKP requires a numerical evolution
of the orbit described by a coupled system of ordinary differential equations
(ODEs) and a root-finding method to solve the Kepler problem, and therefore
lacks the computational efficiency required for most data analysis
applications. The post-circular (PC) formalism [910] provides a method to
recast the time-domain response function $h(t)$ into a form that permits an
approximate fully analytic Fourier transform in stationary phase approximation
(SPA), under the assumption that the eccentricity is small, leading to non-
spinning, eccentric Fourier-domain inspiral waveforms as a simple extension to
the quasi-circular PN approximant TaylorF2. More recent work has extended the
PC formalism to 3PN, with a bivariate expansion in eccentricity and the PN
parameter [911, 903], and has included previously unmodeled effects of
periastron advance [912].

The parameter space coverage of Fourier-domain waveforms in the PC formalism
is limited by the necessary expansion in small eccentricity. Newer models aim
for validity in the range of moderate to high eccentricities, by utilizing
numerical inversions in SPA and resummations of hypergeometric functions to
solve orbital dynamics [913, 750] or by applying Padé approximation on
analytic PC schemes expanded into high orders in eccentricity. A semi-analytic
frequency-domain model for eccentric inspiral waveforms in the presence of
spin-induced precession has been developed with the help of a shifted uniform
asymptotics (SUA) technique to approximate the Fourier transform, a numerical
treatment of the secular evolution coupled to the orbital-averaged spin-
precession, and relying on small eccentricity expansion only for the amplitude
[914].

#####  4.2.3.4 Insight from scattering

One can gain surprising insight in the relativistic two-body problem by
investigating unbound orbits. Studying such systems implies analyzing the
scattering of compact objects and the approximation to the two-body problem
most naturally applicable to it, the PM expansion. Recently, there has been
renewed interest in the subject, as it has been realized that PM information
from unbound systems can be transferred to bound ones, as done for instance
via Hamiltonians [915, 916, 917] or between gauge-invariant quantities [918,
919, 920, 921, 922]. Moreover, PM expansions can be independently obtained
from scattering-amplitude calculations [923], as done at 3PM order [924, 925,
926, 927] for the nonspinning sector (see also earlier results at 2PM order
[928, 929, 930, 931, 932] and at 4PM [933, 934, 935, 936, 937, 938]). Another
method that has been successfully used is the worldline quantum field theory
approach [939, 940, 941].

The 3PM radiative contribution, which cancels a divergence in the 3PM
conservative part, has been obtained in Refs. [942, 943, 944, 945, 946, 947,
948]. The radiative effects in PM expansions have been further explored in
[949, 950, 951, 952] and the radiative contributions to the scattering problem
were obtained in a PN expansion [953, 954, 955, 956, 957, 958, 959].

A Hamiltonian including spin-orbit coupling at the 2PM approximation has been
obtained [960] by extracting the dynamical information from a scattering
situation with the help of “scattering holonomy” [961, 919]. Likewise effects
at ‘spin-squared’ have been determined at the 2PM level using both amplitude
[962] and the EFT formalism [963, 964, 965, 966, 967, 968]. To 3PM order, SO
and SS contributions have been determined [969]. Effects of arbitrary orders
in spin in the scattering angle have also been studied by means of purely
classical methods [970, 877, 971]. Those results have been recovered at 1PM
order through the computation of quantum scattering amplitudes for minimally
coupled massive spin-$n$ particles and gravitons in the classical limit, as
$n$ goes to infinity [972]. Notably, the tree-level amplitude has been found
to generate the full series of black-hole spin-induced multipole moments
[973]. Tidal effects in PM expansions have also been investigated [974, 975,
976, 977, 978, 979, 980, 981].

The aforementioned results have naturally prompted discussions and comparisons
between the GW astrophysics and scattering-amplitude communities, which is
made possible by the use of pivotal gauge-invariant quantities for
comparisons. The scattering angle of unbound compact-object interactions is
one such pivotal quantity that, at least in a perturbative sense through the
orders so-far considered, is thought to encapsulate the complete conservative
dynamics. Not only does it provide a common ground to exchange information
between independent PM calculations [982] or between PM and EOB schemes [915,
916, 970, 983, 984, 985, 986, 987], but it has also proven extremely useful to
PN theory.

The PM expansion of the scattering angle allows one to extract previously
unknown PN information through a synergistic combination of constraints from
its simple mass-ratio dependence and gravitational self-force (GSF) results,
as first realized at 5PN [988, 989]. Such construct has been exploited to
partially calculate the 6PN dynamics [818, 990], allowing a nontrivial check
at 3PM-6PN order of the 3PM result [925, 991] (see also [992, 927]), as well
as to the generic spin-orbit and aligned bilinear-in-spin sectors at NNNL
[993, 994] (see also [995] for more on the interface between PM and GSF
theory). GSF scattering observables can also provide a powerful handle on PM
dynamics across all mass ratios. Calculations of the scattering angle to
first-order (second-order) in the mass ratio fully determine the complete two-
body Hamiltonian through 4PM (6PM) order [996].

#####  4.2.3.5 First laws and gauge invariant comparisons.

The orbital dynamics of a binary system of compact objects exhibits a
fundamental property, known as the first law of binary mechanics, that takes
the form of a simple variational relation. This formula relates local
properties of the individual bodies (e.g. their masses, spins, redshifts, spin
precession frequencies), to global properties of the binary system, (e.g.
gravitational binding energy, total angular momentum, radial action variable,
fundamental frequencies). The first law was first established for binary
systems of nonspinning compact objects moving along circular orbits [997], as
a particular case of a more general variational relation, valid for systems of
black holes and extended matter sources [998]. This first law was later
extended to generic bound eccentric orbits [999], including the effect of the
GW tails that appear at the leading 4PN order [1000], as well as to spinning
compact binaries, for spin aligned or anti-aligned with the orbital angular
momentum [1001]. In the context of the small mass-ratio approximation (see
Sec. 4.3), analogous relations were established for a test particle or a small
self-gravitating body orbiting a Kerr black hole, by accounting for the
conservative part of the first-order gravitational self-force (GSF) [1002,
1003]. These various first laws of binary mechanics have proven useful for a
broad variety of applications, including to:

  * $\bullet$

Determine the numerical values of the “ambiguity parameters” that appeared in
the derivations of the 4PN two-body equations of motion [1004, 802, 803, 804,
806, 807, 808, 809];

  * $\bullet$

Compute the exact first-order conservative GSF contributions to the
gravitational binding energy and angular momentum for circular-orbit
nonspinning black hole binaries, allowing for a coordinate-invariant
comparison to NR results [1005];

  * $\bullet$

Calculate the GSF-induced correction to the frequency of the Schwarzschild
[1005, 1006] and Kerr ISCO [1007, 1008];

  * $\bullet$

Calibrate the potentials that enter the EOB model for circular orbits [1009,
1006] and mildly eccentric orbits [1010, 1011, 1012], and spin-orbit couplings
for spinning binaries [1013];

  * $\bullet$

Test the cosmic censorship conjecture in a scenario where a massive particle
subject to first-order GSF falls into a nonspinning black hole along unbound
orbits [1014];

  * $\bullet$

Define the analogue of the redshift variable of a particle for black holes in
NR simulations, this allowing further comparisons to the predictions of the PN
and GSF approximations [1015, 1016];

  * $\bullet$

Provide a benchmark for the calculations of the first-order GSF-induced
frequency shift of the Schwarzschild innermost bound stable orbit [1017] and
the second-order GSF contribution to the gravitational binding energy [1018].

The first law of binary mechanics for circular orbits has been extended to
account for finite-size effects such as the rotationally-induced and tidally-
induced quadrupole moments of the compact objects [1019, 1020]. Other
interesting directions would be to extend the first law to generic precessing
spinning compact binaries, and to establish it in the context of the PM
approximation and for unbound orbits.

####  4.2.4 Environmental effects

The evolution of relativistic magneto-hydrodynamics, circumbinary disks around
supermassive binary black holes have been studied in the weak-field
approximation. This relies on the techniques of matched asymptotics to
analytically describe a non-spinning binary system up to 2.5PN order [1021],
on which the magnetic field is then numerically evolved [1022, 1023]. These
results have later been extended to spinning, non-precessing binary black-
holes [1024, 1025] and to spinning, precessing black-hole spacetimes [1026].

Another astrophysical environmental effect concerns the presence of a third
body around a massive binary back hole. Hierarchical triplets may undergo
several type of resonances [1027], that can result for example in Kozai-Lidov
oscillations [1028, 1029]. This is an interchange between the eccentricity of
the two-body inner orbit and its inclination relative to the plane of the
third body studied. Such an effect is particularly interesting for LISA as it
involves high eccentricity systems and can result in an enhancement of
gravitational radiation. A general-relativistic treatment of triple
hierarchical systems has been performed in a weak-field approximation up to
2.5PN order in the dynamics [1030]. The quadrupolar and octupolar waveforms
have also been obtained by numerically integrating the three-body system
trajectories from the PN expansion [1031]. Other studies on the Kozai-Lidov
oscillations in GR have employed an expansion in powers of the ratio between
the two semi-major axis up the the hexadecapole order [1031] and later to
second-order in the quadrupolar perturbation [1032]. Finally, an effective
two-body approach to the hierarchical three-body problem has been proposed
[1033] and the mass quadrupole to 1PN order was derived under this formalism
[1034].

####  4.2.5 Challenges

Post-Newtonian information in data analysis enters primarily through the use
of semi-analytical models (see Secs. 4.5,4.6). To avoid parameter estimation
biases due to waveform uncertainties, developments in PN theory will be needed
to increase the accuracies of these models for the next generation of GW
detectors (see 3 and [1] for a discussion in the context of third-generation
ground-based detectors).

One of the most significant challenges in the PN formalism will be the
completion of the 5PN order. First, due to the so-called “effacement
principle” [1035], the point-particle approximation breaks at this order.
Finite-size and tidal effects have to be taken into account, the latter being
known up to 7PN [1036, 1037, 1038, 1039]. In addition, due to the complexity
of the computations involved, it seems highly unlikely to achieve the
derivation of the complete gravitational phase at 5PN during the next decade,
even if partial results are already obtained. In this regard, the interplay
with GSF and scattering amplitudes will be crucial to make significant
progress in this direction.

Exploiting the links between asymptotic symmetries and hereditary effects in
the PN formalism should improve the matching to NR as memory effects are
starting to be included in NR waveforms [682] (see also Sec. 4.1).

Regarding spinning binary systems, although the orbital dynamics are known up
to the 4PN order, it might not be sufficient for LISA data analysis. In the
near future one should primarily focus, at lower orders, on the exterior
gravitational field, outside the matter source, extending to infinity. Indeed,
new multipole moments will have to be computed to reach the same level of
accuracy for the GW amplitude as for the flux, notably the 3PN SS
contributions to the current quadrupole. The next step will consist in moving
on to 4PN order, to do as well as the current equations of evolution. For this
purpose, the NNL SS piece of the 4PN mass quadrupole will be required.

Another important topic, in order to combine information coming from the PN
framework and numerical relativity, will be to connect the magnitude of the
spin used in PN expressions to the rigorous definitions employed in numerical
simulations [1040, 1041]. The problem of comparing the spin axis is even more
intricate since, for now, its direction has not been given a satisfactory
unambiguous meaning in any of these schemes [1042].

Regarding eccentric waveforms, most current models have focused on a low
eccentricity expansion and non-spinning or aligned spin systems. It is crucial
to make progress in both directions. First one has to go beyond the small
eccentricity approximation, as some have started investigating [913, 750].
Second, it will be important to have reliable fully precessing and eccentric
waveforms in order to span all the parameter space expected for LISA sources.

Notwithstanding the amount of work and progress already made, studies of the
scattering of compact objects still explore a relatively uncharted territory.
Obvious extensions of the work reported above will involve pushing scattering
calculations of both spinning and nonspinning systems to higher PM orders in
the conservative sector, both from a fully classical and amplitude approach.
Such new information would help in obtaining new nontrivial PN information via
the continued exploitation of the scattering angle’s mass-ratio dependence, as
well as continuing the fruitful exchange between the communities. The
interconnectedness of the PM, PN and GSF approximations should also be
explored in the context of tidal effects (information from PM is available,
e.g., [974]). Scattering studies should also explore dissipative radiation-
reaction effects. It is important to push these calculations, so that the
accuracy gained in the dissipative dynamics goes on par with that of the
conservative sector.

###  4.3 Small-mass-ratio approximation (gravitational self-force)

Coordinators: Marta Colleoni and Adam Pound
Contributors: S. Akcay, E. Barausse, B. Bonga, R. Brito, M. Casals, G.
Compere, A. Druart, L. Durkan, A. Heffernan, T. Hinderer, S. A. Hughes, S.
Isoyama, C. Kavanagh, L. Kuchler, A. Le Tiec, B. Leather, G. Lukes-
Gerakopoulos, O. Long, P. Lynch, C. Markakis, A. Maselli, J. Mathews, C.
O’Toole, Z. Sam, A. Spiers, S. D. Upton, M. van de Meent, N. Warburton, V.
Witzany

####  4.3.1 Description

When the secondary object in a binary is significantly smaller than the
primary, we can treat the mass ratio $\epsilon=m_{2}/m_{1}=1/q$ as a small
parameter and seek a perturbative solution for the spacetime metric,
$g_{\alpha\beta}=g^{(0)}_{\alpha\beta}+\epsilon
h^{(1)}_{\alpha\beta}+\epsilon^{2}h^{(2)}_{\alpha\beta}+\ldots$ The background
metric $g^{(0)}_{\alpha\beta}$ then describes the spacetime of the primary in
isolation, typically taken to be a Kerr BH. At zeroth order, the secondary
behaves as a test mass in the background, moving on a geodesic of
$g^{(0)}_{\alpha\beta}$. At subleading orders, it generates the perturbations
$h^{(n)}_{\alpha\beta}$, which then exert a self-force back on it,
accelerating it away from geodesic motion and driving the inspiral. These
perturbations $h^{(n)}_{\alpha\beta}$ also encode the emitted waveforms.

The perturbative field equations and trajectory of the secondary are
determined using GSF theory, a collection of techniques for incorporating a
small gravitating object into an external field [1043, 1044, 1045, 98, 1046].
Using matched asymptotic expansions (or EFT [1047]), the small object is
reduced to a point particle, or more generally to a puncture in the spacetime
geometry, equipped with the object’s multipole moments. For a compact object,
higher moments scale with higher powers of the small mass $m_{2}$, such that
one additional multipole moment appears at each order in $\epsilon$. This
multipolar particle (or puncture) is found to obey a generalized equivalence
principle, behaving as a test body in a certain effective, smooth, vacuum
metric $g^{(0)}_{\alpha\beta}+\epsilon
h^{{\mathrm{R}}(1)}_{\alpha\beta}+\epsilon^{2}h^{{\mathrm{R}}(2)}_{\alpha\beta}+\ldots$
[1048, 1049, 1050, 1051, 1052, 1053, 1054, 1055, 1056, 1057, 1058, 1044, 1059,
1060].

The regular fields $h^{{\mathrm{R}}(n)}_{\alpha\beta}$ that influence the
secondary’s motion can be calculated directly, using a puncture scheme in
which they are the numerical variables [1061, 1062]. Alternatively, they can
be extracted from the full fields $h^{(n)}_{\alpha\beta}$ using a mode-by-mode
subtraction [1063, 1064, 1065, 1066, 1067, 1068, 1069] or other methods [1070,
1071], as reviewed in [1072, 1073, 1046]. In either approach, the perturbative
field equations can be solved using the methods of BHPT [1074, 1075, 1076,
1077, 1078, 1079, 1080, 1081, 1082, 1083, 740, 1046].

Figure 5:  Geodesics of Kerr spacetime. Left: a non-resonant geodesic, which
fills a toroidal shape. Right: a resonant geodesic. Images taken from [98].

A key feature of this model is its clean separation of time scales: because
the self-force is small, inspirals occur slowly, over $\sim 1/\epsilon$
orbits. On the time scale of a few orbits, the secondary’s trajectory is
approximately a bound geodesic of $g^{(0)}_{\alpha\beta}$ (illustrated in Fig.
5), which in Kerr spacetime is generically triperiodic, undergoing radial,
polar, and azimuthal motion with frequencies $\Omega_{r}$, $\Omega_{\theta}$,
and $\Omega_{\phi}$ [1084, 1085, 1086, 1087, 1046]. Over the long inspiral,
the frequencies slowly evolve due to dissipation. The field equations can
therefore be solved using a two-timescale expansion [1088, 1089, 1090, 1091,
1092, 1093, 1046],
$h^{(n)}_{\alpha\beta}=\sum_{k^{A}\in\mathbb{Z}}h^{(n,k^{A})}_{\alpha\beta}(J^{B},x^{a})e^{ik^{A}\psi_{A}}$.
Here $x^{a}$ are spatial coordinates, and all time dependence is encoded in
the set of system parameters $J^{A}$ (the secondary’s orbital energy, angular
momentum, and Carter constant, the primary’s mass and spin, etc.), and the set
of orbital phases $\psi_{A}$. The parameters $J^{A}$ evolve slowly, on the
inspiral time scale $t\sim 2\pi/(\epsilon\Omega)$, while the phases $\psi_{A}$
evolve on the orbital time scale $t\sim 2\pi/\Omega$, with evolution equations
of the form [1089, 747, 1092, 1093, 1046]

$\displaystyle\frac{dJ^{A}}{dt}=\epsilon
G_{(1)}^{A}(J^{B})+\epsilon^{2}G^{A}_{(2)}(J^{B})+O(\epsilon^{3}),\quad\frac{d\psi_{A}}{dt}=\Omega^{(0)}_{A}(J^{B})+\epsilon\Omega^{(1)}_{A}(J^{B})+O(\epsilon^{2}).$
(13)

Formulating the problem in this way enables practical methods of solving the
Einstein equations [1092, 1046] as well as facilitating rapid waveform
generation [747, 1046, 1094].

The two-timescale description will frame much of the discussion below.
However, there are some special cases of interest that it does not apply to,
such as scattering orbits [1095, 1096, 1097, 1098, 1099]. Even in an inspiral,
there are regions of parameter space in which the approximation breaks down:
at the end of the inspiral, when the secondary transitions into a final plunge
into the primary [1100, 1101, 1102, 1103, 1104, 1105], and during resonances,
which occur when at least two of the frequencies $\Omega_{A}$ have a rational
ratio, causing a linear combination of phases, $k^{A}\psi_{A}$, to become
approximately stationary [1106, 1107, 1108]. Such resonances can arise from a
variety of physical causes, outlined below, with significant observational
consequences.

Many of the core tools in GSF modelling, as well as advanced codes described
in the sections below, have been consolidated in the open-source Black Hole
Perturbation Toolkit [740], which provides a hub for GSF code development, as
well as in the Black Hole Perturbation Club [1109]. In particular, the Fast
EMRI Waveforms (FEW) package [1110], which exploits GSF models’ multiscale
structure, has provided a flexible framework for rapid waveform generation.

####  4.3.2 Suitable for what sources?

Historically, the GSF approximation has been motivated by EMRIs, with mass
ratios in the interval $10^{-7}\lesssim\epsilon\lesssim 10^{-3}$ as far as
LISA sources are concerned (see Table 3). The GSF model’s accuracy for a
typical EMRI can be estimated from the evolution equations (13). On the
inspiral time scale, the solutions for the orbital phases, and therefore the
GW phase, take the form

$\psi_{A}=\frac{1}{\epsilon}\left[\psi^{(0)}_{A}+\epsilon\psi^{(1)}_{A}+O(\epsilon^{2})\right]\\!.$
(14)

Following Ref. [1089], the leading term in this expansion is referred to as
adiabatic order (0PA), and the $n$th subleading term as $n$th post-adiabatic
order ($n$PA). An adiabatic approximation, which has large phase errors
$\sim\epsilon^{0}$ over an inspiral, is expected to suffice for detection of
EMRIs using the semi-coherent searches discussed in Sec. 3.3.1. A 1PA
approximation, which will have phase errors $\sim\epsilon\ll 1$ rad over the
final year of inspiral, should suffice for parameter extraction even for loud
EMRIs with SNR $>50$.

A large and growing body of evidence [1111, 1112, 1113, 1114, 1005, 690, 1115,
1116, 1008, 1015, 1016, 1117, 506, 1118, 1119, 505, 717] suggests that the GSF
approximation can also accurately model IMRIs (see Table 4). It can even
provide a useful model of CO binaries with comparable masses when it is re-
expressed as an expansion in powers of the symmetric mass ratio $\nu\equiv
m_{1}m_{2}/(m_{1}+m_{2})^{2}=\epsilon+O(\epsilon^{2})$. We refer to Ref.
[1120] for an early review and Refs. [505, 717] for recent analyses of the
accuracy of 1PA models in the $1\lesssim q\lesssim 100$ regime. Further study
is required to assess whether 1PA models meet expected accuracy requirements
for specific classes of astrophysical sources in this range of mass ratios.

Besides compact binary inspirals, the GSF approximation can be used for other
classes of sources. A leading-order GSF model, comprising a point source on a
geodesic trajectory, should be suitable for XMRIs [72, 1121] and fly-by burst
signals [1095, 1096, 1097, 1098, 1099]. GSF theory may also be relevant for
modelling gravitational waves from cosmic strings [285, 286, 287]. Inspirals
of less compact bodies into MBHs (e.g., white or brown dwarfs or main-sequence
stars prior to tidal disruption), and three-body “binary EMRIs”, can be
modelled by including sufficiently high multipole moments in the particle or
puncture [1122, 87].

####  4.3.3 Status

The current status of GSF models is summarized in Fig. 6, which also
summarizes the necessary inputs for a waveform model at 0PA and 1PA order.

Figure 6:  Progress in modelling EMRIs using GSF methods. ‘1SF’ and ‘2SF’
indicate calculations involving the first or second-order self-force,
respectively, and ‘spin effects’ indicates calculations that take into account
the secondary’s spin (the spin of the primary is accounted for in all ‘Kerr’
calculations. ‘Snapshot calculations’ (single tick) are ones in which self-
force effects are calculated on fixed geodesic orbits.

#####  4.3.3.1 Adiabatic approximation

At leading order in the two-timescale expansion, only the geodesic frequencies
$\Omega^{(0)}_{A}$ and time-averaged dissipative piece $G_{(1)}^{A}$ of the
first-order self-force are required to drive the evolution [1085, 1123]. The
secondary’s motion in this approximation can be regarded as an adiabatic
inspiral through a sequence of geodesic orbits [1089, 1124]. The waveform can
then be built from a corresponding sequence of “snapshot” Fourier mode
amplitudes together with the leading phases in Eq. (14). Alternatively,
outside the two-timescale scheme, a time-domain field equation can be solved
with an adiabatically inspiraling source particle [1125, 1126, 1127, 1128,
1117].

$G_{(1)}^{A}$ is most conveniently computed using “flux-balance” formulas
[1129, 1130, 1131, 1132, 1133], in which the “work” done by the dissipative
self-force balances the fluxes of gravitational waves to infinity and down the
horizon (or flux-like quantities in the case of the Carter constant). Methods
of computing fluxes in GSF theory are well developed [1134, 348, 1135, 1136,
1137, 1138, 1125, 1139, 1128], and have been numerically implemented for
generic (inclined and eccentric) geodesics about a Kerr BH [1140, 1141, 1142,
1143]. The fluxes have also been calculated analytically (by expanding
$h^{(1)}_{\alpha\beta}$ in a PN series [1144, 1145, 1081]) to high PN order
for circular [1146, 1147, 1148], eccentric [1149, 1150, 1151], and generic
orbital configurations [1152, 1153, 1154].

After decades of progress [1155, 70, 1156, 1157, 1158, 1159, 1103, 1160,
1161], adiabatic inspirals and waveforms are now being computed for generic
orbits in Kerr, both numerically [1162] and within the analytical BHPT-PN
framework mentioned above [1163].

Once fluxes and waveform mode amplitudes have been computed across the
parameter space, adiabatic waveforms can be efficiently generated using a
combination of neural network and reduced-order techniques [1160, 1110], a
method that lends itself to graphical processing unit (GPU) acceleration
techniques; see Sec. 5.2. This approach should ultimately meet the efficiency
requirements described in Sec. 3.3.1. However, due to the high-dimensional
parameter space, significant work remains to populate the parameter space. As
a consequence, data analysis development has so far relied on semi-
relativistic “kludge” models that approximate the adiabatic evolution and
waveform generation [1164, 1165, 1166, 488, 1167].

#####  4.3.3.2 Post-1-adiabatic self-force effects: 1SF

The 1PA terms in the two-timescale expansion ($\Omega^{(1)}_{A}$ and
$G^{A}_{(2)}$) take as input the conservative piece of the first-order self-
force (1SF, $\propto h^{(1)}_{\alpha\beta}$) and the dissipative piece of the
second-order self-force (2SF, $\propto h^{(2)}_{\alpha\beta}$). Once these
ingredients have been computed, a 1PA waveform-generation scheme takes the
same form as a 0PA one [1092, 1046].

Numerical computations of the full first-order self-force have made tremendous
progress, evolving from Lorenz-gauge calculations [1083, 1168, 1169, 1170,
1171, 1172, 1007, 1073] to more efficient methods relying on radiation gauges
[1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180]. These developments
culminated in 1SF calculations along generic bound geodesics in Kerr spacetime
[1181], and work is still ongoing to further improve 1SF methods [1182, 1183,
1184].

High-precision numerics and BHPT-PN methods have also made possible high-order
PN expansions of numerous conservative invariants, such as the Detweiler
redshift [1185, 1186, 1187, 1188, 1177, 1189, 1190, 1191], the periastron
advance [1008, 1192], and other invariants [1193, 1194, 1195, 1196, 1197,
1198], which have played a key role in the synergies with PN theory and EOB
described in Secs. 4.2 and 4.5.

Modern 1SF calculations mostly rely on frequency-domain methods and mode-by-
mode subtraction or puncture schemes to calculate
$h^{\mathrm{R}(1)}_{\alpha\beta}$. However, there have also been advances in
time domain (TD) calculations based on either finite-difference schemes [1199,
1124, 1200, 1097, 1201] or spectral methods [1202, 1203, 1204, 1205, 1206,
1207, 1201, 1208], along with improved time-stepping methods [1209, 1210,
1201, 1211]. There is also ongoing work to directly obtain the retarded Green
function [1212, 1071, 1213, 1214, 1215], which would then allow direct
evaluation of the self-force [1048, 1049].

As at 0PA, calculations of the self-force and of GW amplitudes across the
parameter space have been used to simulate self-forced inspirals [1216, 1217,
747, 1094, 1218] and generate waveforms [1219, 1217]. To date, this has only
been done for equatorial orbits and quasi-spherical orbits [1218], due to the
computational expense of current 1SF calculations for generic orbits.

#####  4.3.3.3 Post-1-adiabatic self-force effects: 2SF

The dissipative piece of the second-order self-force ($G_{(2)}^{A}$)
contributes to the GW phasing at the same 1PA order as the first-order
conservative self-force, but calculations of it are less mature. After years
of development of the governing formalism [1220, 1056, 1057, 1058, 1221, 1222,
1059, 1060] and of practical implementation methods [1223, 1224, 1073, 1091,
1225, 1226, 1046, 1227, 1228], Refs. [1018, 1118, 1119] recently carried out
the first concrete calculations of physical second-order quantitities in the
restricted case of quasicircular orbits around a Schwarzschild BH. These
calculations culminated in the first complete 1PA waveforms in Ref. [1119].
Figure 7 shows a comparison between one of these 1PA waveforms and an NR
waveform for $q=10$. As alluded to in Sec. 4.3.2, the 1PA approximation agrees
well with NR even at this moderate mass ratio.

These 1PA waveforms have also more recently been extended to include a slowly
spinning primary [1229]. Sec. 4.3.5 discusses the next major barrier in GSF
calculations: 2SF calculations with a rapidly spinning primary, eccentricity,
and inclination.

Figure 7: 1PA GSF waveform for a quasicircular, nonspinning binary with mass
ratio $q=10$ (orange). The inset shows a zoomed portion of the waveform near
the merger. Also included for comparison are the 0PA GSF waveform (blue, inset
only) and the waveform for the same binary produced using an NR simulation in
the SXS catalog (SXS:BBH:1107, in black). The three waveforms are aligned in
time and phase at $t=320M$, when the orbital separation is $\approx 13.83M$.
Image reproduced from Ref. [1119].

#####  4.3.3.4 Transient resonances

Resonances are a ubiquitous feature of the strong-field dynamics around BHs in
GR [1230, 1231], and their observational imprints on waveforms can be
significant. Transient self-force resonances, between the orbital frequencies
$\Omega_{r}$ and $\Omega_{\theta}$, will have strong observational
consequences for EMRIs [1123, 1106, 1232, 1233, 1107]. Essentially all LISA-
type EMRIs pass through at least one dynamically significant resonance,
leading to a large, $O(1/\sqrt{\epsilon})$ (i.e., 0.5PA) contribution to the
waveform phase [1234, 71, 1235]. The magnitude of this effect depends
sensitively on the orbital phase at the resonance. As a result, modelling it
requires 1PA accuracy prior to resonance. There are ongoing efforts to
understand the impact of these resonances and include them in evolutions
[1236, 1237, 1238, 1132, 1239, 1133, 1163, 1108, 1240, 1241, 94].

There also occur resonances between the $r$ and $\varphi$ frequencies and
between the $\theta$ and $\varphi$ frequencies. These do not change the
intrinsic inspiral dynamics, but they can lead to a strong net emission of
linear momentum that results in a ‘kick’ to the system’s center of mass [1242,
1243]. The maximum kick velocity can reach $\sim 30,000\times\epsilon^{3/2}$
km/s (depending on the primary spin and on the orbital eccentricity), which,
for IMRI systems, could be comparable to the escape velocity from a typical
globular cluster.

#####  4.3.3.5 Merger and ringdown

The two-timescale expansion breaks down as the secondary object approaches the
separatrix between bound and plunging orbits. There, the secondary enters a
gradual transition across the separatrix [1100, 1101, 1103, 1104] followed by
an approximately geodesic plunge. The transition motion has been mostly
studied within the EOB framework using resummed PN expansions [1100, 1244,
1245, 1246, 1247, 1248, 1249, 1250, 1251] or within the GSF approach using
simplifying approximations [1101, 1252, 1253, 1127, 1102, 1104, 1103]. Merger-
ringdown waveforms have been generated from the final plunge at first order
[1254, 1255], and complete inspiral-merger-ringdown waveforms have been
generated by solving time-domain equations for first-order perturbations
sourced by the combined inspiral-transition-plunge motion [1126, 1127, 1128,
1117, 766]. Only recently, Refs. [1256, 1105, 1257] took a first step toward
1PA/2SF inspiral-merger-ringdown waveforms by developing a systematic
expansion of the transition motion and metric perturbation that matches to the
two-timescale expansion for quasicircular, equatorial inspirals into a Kerr
BH.

#####  4.3.3.6 Spin and finite-size effects

When its multipole structure is accounted for, the secondary obeys the
Mathisson-Papapetrou-Dixon equations of a test body [841, 843, 1258, 1259,
844] in the effective metric; such a result is expected to hold to all
perturbative orders for both BHs and material bodies [1260, 1055, 1044].

The spin contributes first-order conservative and second-order dissipative
forces in the equations of motion, leading to 1PA, $\mathcal{O}(\epsilon^{0})$
contributions to the final inspiral phase [1261]. Spin also generally breaks
the integrability of Kerr geodesic motion [1262]; the resulting potential
phenomena, like prolonged resonances, can leave a significant imprint on the
gravitational waves [1263, 1264, 1265, 85, 1266] and may contribute
$\mathcal{O}(1)$ to the GW phase. Away from these resonances, the spinning-
particle motion is perturbatively separable to linear order in spin [85] and
easily incorporated into a two-timescale approximation [1267, 1261].

For dissipative effects, a spin-flux balance law has been established [1268],
and GW fluxes were computed for a variety of binary configurations [1269,
1270, 1271, 1272, 1273, 1274, 1275, 1268, 1266, 1276, 1277]. Waveforms from
generic inspirals into a Schwarzschild BH have also been computed including
1SF and first-order conservative spin effects [1278] but excluding dissipative
spin effects. Waveforms including all 1PA spin effects from circular
equatorial inspirals into Schwarzschild [1261] and Kerr BHs [131] have been
produced, and the spin’s complete 1PA contribution to the GW phase has been
calculated for eccentric equatorial inspirals into a Kerr BH [1279]. A
formulation of the first-order conservative spin-forced motion for generic
orbits about a Kerr BH [1280, 1281] paved the way for calculations of the
spin’s complete impact on fully generic 1PA waveforms, beginning with recent
snapshot computations of energy and angular momentum fluxes [86].

The secondary’s higher moments are unlikely to contribute at 1PA order. The
tidal quadrupole’s contribution to the acceleration scales with $\epsilon^{4}$
[1282] and can only become relevant for non-compact objects [1267]. The spin
also induces a quadrupole deformation of the secondary, which creates an
acceleration $\sim\epsilon^{2}$ [1122, 87] that is unlikely important for
EMRIs (up to resonances) but might be relevant for IMRIs [87].

####  4.3.4 Environmental effects

In the preceding sections, we have focused on the simplest GSF model: an
isolated vacuum BH orbited by a single body. This will not fully describe
astrophysical small-mass-ratio binaries. Accreting matter will be present at
some level, and there are likely to be additional bodies nearby. Accounting
for these effects will be necessary to achieve the stringent accuracy
requirements that EMRIs, in particular, place on GSF models.

There are three ways we can modify GSF models to include such effects: (i)
perturbatively, by adding “small” matter fields or metric perturbations to the
spacetime, which then exert small forces on the secondary; (ii) non-
perturbatively, by modifying the background spacetime and introducing “large”
matter fields; (iii) by changing boundary conditions, imagining (for example)
modifications very near the BH horizon or non-asymptotically flat
perturbations due to external matter. Sec. 6 discusses how these classes of
modifications can also be used to incorporate beyond-GR effects.

Although important aspects of the astrophysical environment of EMRIs are
uncertain, it is expected that there will be other stellar-mass bodies nearby.
These bodies create an additional metric perturbation, which should be
sufficiently small to be treated as a linear perturbation. The influence of
the perturbation is generally negligible, except at moments when two or more
of the orbital frequencies characterising the perturbed EMRI system become
commensurate. At those times, an EMRI experiences a tidal resonance [1283], an
orbital resonance akin to the mean motion resonance known from planetary
dynamics [1284]. Calculations across the parameter space showed that a single
resonance can dephase the waveform by several radians over the inspiral [89,
1285, 94, 1286], and that most EMRIs will cross multiple tidal resonances
before plunge. If LISA can reliably measure these resonances, they can be used
to learn about the tidal environment of EMRI systems.

Accretion and gas interactions will also be present at some level in any
astrophysical EMRI system. Order-of-magnitude estimates suggest they will be
negligible unless the primary BH powers an AGN (in which case it would be
surrounded by an accretion disk with which the secondary is likely to
interact) [64, 93, 1287, 1288]. Since AGNs are believed to make up 1–10% of
local galaxies, only a comparable fraction of EMRIs are expected to be
significantly affected by accretion and gas interactions. For those systems,
accretion onto the primary and secondary, as well as dynamical friction from
the disk and planetary-like migration within it, are expected to cause secular
effects comparable to those of GW fluxes [64, 93, 1287]. Recently work has
begun to incorporate the disk and other potential nonvacuum effects into GSF
models by adding torques to 0PA Kerr models [1289] or by working
perturbatively on an exact nonvacuum background [1290, 1291].

The direct gravitational pull from AGN disks is likely to be negligible [64],
but if more dense disks/rings exist in nature, the effect may be significant
[1292]. Exact spacetime solutions describing thin disks or rings around BHs
[1293, 1294, 1295, 1290] show that orbits in such situations are non-
integrable, exhibiting characteristic phenomena like chaos and prolonged
resonances [1296]. If these effects are significant in realistic scenarios,
GSF models should be able to account for them through perturbative corrections
on top of a Kerr background [1297, 1298, 94, 1299].

####  4.3.5 Challenges

The principal goal of GSF waveform modelling is to develop complete 1PA models
for generic orbital configurations around a spinning, Kerr BH. To be complete,
these models must include the spin of the companion, transitions across
resonances, and the final plunge. They must also be sufficiently modular to
incorporate beyond-GR and environmental effects, and there are strong
motivations to explore other regions of the parameter space, such as scatter
orbits and comparable masses. We summarize here the main challenges in
developing and implementing such models.

#####  4.3.5.1 Post-1-adiabatic calculations

1PA models are currently missing two ingredients: the effects of the
companion’s spin for generic binary configurations, and dissipative 2SF
effects; see Fig. 6. Both must be incorporated into a unified two-timescale
expansion of the field equations, the orbital motion, and the spin evolution.

While substantial work remains to calculate 1PA spin effects for generic
binary configurations, these calculations can leverage existing methods for
point-particle sources. Therefore, 2SF calculations represent the overriding
obstacle to 1PA accuracy. Practical 2SF calculations have, thus far, been
restricted to quasicircular orbits in a Schwarzschild background [1018, 1118,
1119]. As discussed in Sec. 2.2, the majority of EMRIs are expected to have
significant eccentricity and may have highly precessing orbital planes when
they enter the LISA band, meaning 2SF techniques need to be extended to cover
generic orbital configurations. They must also be extended to the realistic
case of a Kerr background.

Eccentricity and inclination bring multiple challenges. Most recent self-force
calculations have utilized decompositions into Fourier modes (the natural
setting of the two-timescale expansion) and angular harmonics. Eccentric,
inclined orbits can require $\sim 10^{5}$ modes, all of which will couple to
one another in the second-order source. The sum of Fourier modes can also
suffer from poor convergence: in the existing 2SF puncture scheme, the source
for $h^{\mathrm{R}(2)}_{\alpha\beta}$ has finite differentiability on the
puncture’s worldline, leading to slow power-law convergence (the Gibbs
phenomenon). At first order, similar problems were overcome using methods of
extended solutions [1300, 1301, 1302, 1171, 1172, 1177, 1181] that restore
exponential convergence. These have inspired a new scheme, applicable at
second order, known as the method of extended effective sources [1303], which
has been demonstrated in the case of a scalar-field toy model for eccentric
orbits about a Schwarzschild BH. Work now remains to apply it to gravitational
perturbations, both at first and second order, and to orbits in Kerr
spacetime.

The extension to Kerr spacetime brings more challenges. Second-order
calculations in Schwarzschild have so far relied on directly solving the
perturbative Einstein equations, which are not separable in Kerr spacetime. At
first order, this problem was overcome using radiation-gauge methods, in which
the metric perturbation is reconstructed from a solution to the (fully
separable) Teukolsky equation [1076, 1078, 1080]. One path to 2SF calculations
is to extend that method to second order [1304, 1305, 1182, 1306]. The
standard method of metric reconstruction [1078, 1080] fails beyond linear
order, but recent work found an extension to all orders [1305, 1182]. Recent
work on second-order flux-balance laws [1307] also suggests that 1PA rates of
change of energy and angular momentum (though perhaps not the Carter constant)
can be computed directly from a solution to the second-order Teukolsky
equation, without the need for metric reconstruction.

An additional obstacle is that in radiation-gauge implementations, the first-
order metric perturbation has gauge singularities that extend away from the
particle. There is a rigorous procedure to extract physical 1SF quantities
despite these singularities [1175, 1179, 1180], but the singularities become
ill defined in the second-order field equations [1182]. Several avenues are
being explored to resolve this problem [1182, 1183, 1308, 1309].

There are also challenges common to all these 2SF calculations: at second
order, the field equations have a noncompact source that falls off slowly at
large distances [1091] and is burdensome to compute due to the strong
nonlinear singularity at the particle [1225]. Recent work [1060, 1306] has
shown that both problems might be mitigated by using gauges adapted to the
lightcone structure of the perturbed spacetime. However, additional work will
be required to implement these gauge choices in a practical numerical scheme.

#####  4.3.5.2 Covering the parameter space

Even once numerical implementations of all necessary ingredients are
available, spanning the full EMRI parameter space remains a considerable
challenge at both 0PA and 1PA orders. This is due to the high dimensionality
of the parameter space and the high computational burden of self-force
calculations, particularly at second order. Covering the EMRI parameter space
will likely involve a combination of (i) using analytic results to reduce the
region of the parameter space where high-precision interpolation of numerical
data is required, (ii) better interpolation methods, and (iii) improvements in
computational efficiency of current numerical calculations. We address each of
these three below in the context of 0PA inspirals, 1PA inspirals, and waveform
calculations.

Computing adiabatic inspirals requires interpolating the rates of change
($G_{(1)}^{A}$) of the orbital energy, angular momentum, and Carter constant
across the four-dimensional parameter space of the primary spin and three
orbital elements. For an adiabatic model sufficiently accurate to build 1PA
corrections upon, we need to interpolate $G_{(1)}^{A}$ to better than a
relative accuracy of $\epsilon^{-1}$. The central challenge is then
constructing such a high-accuracy interpolant over the 4D parameter space.
Recent work has demonstrated advantages of Chebyshev interpolation for this
purpose [1094]. The region of the parameter space, and the accuracy of the
interpolation of the numerical data, can also be reduced by constructing
global fits informed by analytic results.

At 1PA order we must also compute the change in the mass of the primary during
the inspiral [1092], and so the parameter space grows to five dimensions.
There are also corrections due the spin of the secondary, but fortunately,
these can be added on separately. At 1PA order the accuracy requirements of
the contributions are much lower, at $\sim 10^{-2}-10^{-3}$ relative [1217].
This should allow analytic results to assist in reducing the parameter space
that numerical results need to cover. Such analytical results could be
obtained by extending BHPT-PN calculations. Alternatively, they could be
obtained from EOB dynamics, following the programme in Refs. [1310, 1311,
1312, 1313, 1314, 1315, 1316]. Regions of high eccentricity may also be more
easily covered using advanced time-domain codes [1317, 1205, 1206, 1207, 1199,
1097, 1209, 1201, 1208].

#####  4.3.5.3 Extending the parameter space

The challenges above are centered on “vanilla” regions of the parameter space:
the inspiral phase away from resonances, which is amenable to a two-timescale
expansion. However, it is also critical to include accurate transitions across
resonances, particularly in the EMRI regime. In addition to the dominant
$r$-$\theta$ resonances that occur in a self-forced inspiral, prolonged
resonances may occur due to spin, external matter, or a non-Kerr central
object [1233, 1318, 1108]. Similarly, the transition from inspiral to plunge
can also be important, particularly for more comparable masses [1117, 506,
505]. Both resonances and the plunge will require matching the two-timescale
expansion to specialized approximations in those parameter regions [1238, 71,
1046, 1108, 1256, 1105, 94].

There are also reasons to compute self-force effects on scattering orbits,
whether to model hyperbolic, burst sources or to inform PM and PN dynamics, as
described in 4.2.3.4 (and possibly to infer properties of bound, self-forced
orbits [919, 1319]). Explorations of scatter orbits have only recently begun
[1098, 1099, 1320, 1321]. While bound-orbit self-force calculations can
utilise the orbit’s discrete Fourier spectrum, scatter orbits have a
continuous spectrum [1096, 1095], suggesting some clear advantages to
simulations in the time domain [1017, 1098, 1099].

In addition to including more of the two-body GR parameter space, considerable
work must be done to include possible beyond-GR and environmental effects in
GSF models. Fortunately, the GSF model is relatively modular. This means most
additional effects can be added separately, as described in Sec. 6.4, and are
readily incorporated into frameworks such as FEW [1110]. However, modelling
these additional effects will be particularly challenging if they are not
amenable to a two-timescale treatment, or if they are too large to be treated
perturbatively.

###  4.4 Perturbation theory for post-merger waveforms (quasi-normal modes)

Coordinators: Stephen R. Green and Laura Sberna
Contributors: E. Berti, R. P. Macedo, P. Mourier, N. Oshita, M. van de Meent

####  4.4.1 Description

Following a compact binary merger, the remnant object settles into a
stationary state through a process known as the “ringdown”. The ringdown
signal is interesting because it involves a collection of discrete modes that
encode information about the final object. For a BH merger, the no-hair
theorems of general relativity predict that the final state is itself a Kerr
BH (see, e.g., [1322]), yielding precise predictions for the ringdown
frequencies. We consider this case here.

Figure 8:  Real and imaginary parts of the Kerr quasinormal mode frequencies,
coloured by the BH spin (from $a=0$ in dark blue to $a=0.9$ in yellow). Left:
$\ell=m=2$ fundamental and first overtone modes. Right: modes with different
angular numbers $(\ell,m)$.

At late times, the ringdown is best described using BHPT. This involves
expanding the spacetime metric $g_{\alpha\beta}=g^{(0)}_{\alpha\beta}+\epsilon
g^{(1)}_{\alpha\beta}+\epsilon^{2}g^{(2)}_{\alpha\beta}+...$ [1074, 1323,
1075] and solving order-by-order in the Einstein equation. At zeroth order,
the background metric $g^{(0)}_{\alpha\beta}$ describes the remnant Kerr BH.
At first order, $g^{(1)}_{ab}$ satisfies the Einstein equation linearized
about Kerr. Rather than work directly with the metric perturbation, it is
convenient to use instead the Newman-Penrose formalism, expressing the
perturbation in terms of the (complex) Weyl curvature scalars. Equations for
the $s=\pm 2$ linearized Weyl scalars $\psi_{0}$ and $\psi_{4}$ famously
decouple yielding a separable wave-like equation discovered by Teukolsky
[1324, 1076, 1134, 348]. Both $\psi_{0}$ and $\psi_{4}$ uniquely describe
$g^{(1)}_{ab}$ up to gauge and perturbations to other Kerr spacetimes [1077],
and in particular describe the GW degrees of freedom, with $\psi_{4}$ most
relevant for gravitational radiation at infinity. Similarly, spin-weight
$s=\pm 1$ and $s=0$ Teukolsky equations describe electromagnetic and scalar-
field perturbations, respectively.

Homogeneous solutions to the Teukolsky equation can be obtained by imposing
boundary conditions corresponding to the physical requirement that radiation
is infalling at the BH horizon and outgoing at infinity. For each spin weight,
this eigenvalue problem admits a countably infinite, discrete spectrum of
complex-frequency _quasinormal modes (QNMs)_ [1325, 1326], see Fig. 8. Each
QNM is labelled by three integers: two spin-weighted spheroidal-harmonic
quantum numbers used to separate the angular dependence ($\ell\geq 2$,
$-\ell\leq m\leq\ell$) and an “overtone number” $n\geq 0$, which sorts the
frequencies in ascending order of the absolute value of their imaginary part.
Negative (real) frequency modes are related to modes with negative $m$ by the
symmetry $-\omega^{*}_{n\ell m}=\omega_{n\ell-m}$. Kerr quasinormal
frequencies (QNFs) have a strictly negative imaginary part, hence the modes
decay with time, due to absorption by the BH horizon and radiation to
infinity. These frequencies are moreover uniquely determined by the Kerr mass
and spin in agreement with no-hair theorems. For detailed reviews on QNMs,
see, e.g., [1327, 1328, 1329, 1330].

In contrast to normal modes, QNMs of Kerr do not form a complete basis, and
therefore cannot fully describe the ringdown for all times, even at linear
order. At _early_ times, the linear response to a perturbation is dominated by
a “direct” emission of radiation (also known as the prompt response) [1331,
1332, 1328]. At _late_ times, the perturbation is dominated by a power-law
“tail” (associated with a branch cut in the frequency plane of the Green’s
function of the Teukolsky equation), which arises due to the asymptotic
properties of the wave-equation potential at large distances from the horizon
[1332, 1333]. QNMs are most important at _intermediate_ times, and for
realistic inspirals they tend to dominate the majority of the post-merger
signal. They have therefore proven extremely relevant for waveform modeling,
either to analyze the ringdown in isolation or for combined inspiral-merger-
ringdown waveforms. Indeed, a single QNM can be used to infer the mass and
spin of the remnant, whereas additional modes or inspiral information enable
consistency checks on GR [1334, 1335, 612, 1336, 113, 605].

####  4.4.2 Suitable for what sources?

LISA is expected to detect the ringdown of MBHB and IMRI sources for which the
merger occurs during the LISA observational time, and for which remnant QNFs
fall within the LISA frequency band [612, 1337, 1338]. For EMRIs, the ringdown
is likely to be too weak to resolve, a potential exception being those with
near-extremal primaries [1339, 1340].

One can easily estimate the mass range of remnants with ringdowns within the
LISA band. Indeed, the longest lived (fundamental, $\ell=m=2,n=0$) QNM has
frequency $f_{0}\simeq 0.4/(2\pi M_{\mathrm{final}})$, with the exact value
depending on the remnant’s spin. For this to lie within the LISA band,
$[10^{-4},10^{-1}]$ Hz, the mass of the remnant must therefore fall within the
range $[10^{8},10^{5}]~{}\text{M}_{\odot}$. Ringdown modelling will be
particularly important for MBHBs at the higher end of the mass spectrum, since
their inspiral will take place outside the LISA band and contribute little to
the SNR [1341].

####  4.4.3 Status

At linear order in BHPT and for most of the duration of the signal, the
ringdown is well described by a superposition of QNMs. For a given BH mass and
spin, the Kerr QNFs were first computed numerically in [1334]. The most
accurate and reliable numerical method is Leaver’s continued fraction approach
[1342, 1343, 1344] (see Mathematica notebook and data at [1345]) and in
particular its spectral refinement [1346, 1156] (implemented in the Python
package qnm [1347]). Kerr QNFs can also be estimated analytically using WKB
(Wentzel–Kramers–Brillouin) or WKB-inspired approximations [1348, 1349, 1350,
1351, 1352, 1353]. QNM modeling should be performed in a certain preferred BMS
frame [1354, 563]. Afterwards it can be transformed to a different desired
reference frame, e.g., a frame co-precessing with the binary [763, 1355], a
frame set prior to a merger kick [1356], or the post-Newtonian preferred frame
appropriate to modeling the early inspiral [563].

In addition to the frequencies, the QNM component of the ringdown signal also
comprises the amplitude and phase of each mode. Unlike the frequencies, these
quantities are not uniquely determined by the remnant, but rather depend on
the particular initial conditions that led to its formation. In principle, the
amplitudes and phases can be computed from the spacetime right after the
formation of a common BH horizon [1332, 1357], and hence from the binary
parameters [1358, 1359]. However, this data is hard to predict analytically,
and moreover challenging to extract from NR simulations. Common practice is
therefore to calibrate the desired mode amplitudes and phases by fitting
against NR simulations [1360, 1358, 1361, 1362, 1363, 1354, 1364, 1365].
Alternatively, initial data can be estimated at the level of the Teukolsky
equation using an Ori-Thorne procedure [1127, 1366, 1367, 1102, 1368].

For binary BH mergers of comparable masses, the linear ringdown model—a
superposition of a finite number of QNMs with arbitrary amplitudes and
phases—has been extensively compared against NR simulations. At late times
($\gtrsim 10M$ after the peak of the waveform), the general consensus is that
the linear model provides a good, stable, and consistent fit to the numerical
waveform. At these times, the waveform is well described by a small number of
overtones (dominated by the $\ell=m=2$ multipole and overtones with $n\lesssim
2$), and fitting for the frequencies and decay rates allows for accurate
inference of the remnant mass and spin [575, 1362, 1369, 1370, 1371, 1372,
1373, 1354, 1365]. More surprisingly, several studies indicate that the linear
model applies even at much earlier times, provided higher $(\ell,m,n)$ modes
are included in the model [1369, 1370, 1371, 1354]. However, the relevance of
the linear model at these times ($\lesssim 10M$ after the peak of the
waveform) and, in particular, the role of spherical-spheroidal mode mixing,
nonlinear modes and higher overtones ($n>2$), is still actively debated [1374,
1354, 1375, 1376, 1377, 1378, 1379, 1365, 1380, 1381].

A complete linear-order model of the post-merger signal should also include
back-scattering of radiation against the background potential. This is well
approximated by a power-law tail $h\sim t^{-n_{\mathrm{t}ail}}$ at late times,
where $n_{\mathrm{t}ail}=7$ for generic gravitational perturbations of Kerr
[1333, 1382, 1383]. So far the tail contribution has not been confidently
identified in NR simulations of binary BH mergers, so it is usually neglected
in waveform modeling.

At very early times (around the peak amplitude) we expect higher order BHPT to
become relevant. Beyond-linear-order perturbations satisfy the same equations
as at linear order, but now with source terms made up of lower-order
perturbations [1384, 1304, 1385, 1386, 1387, 1388, 1306]. It is therefore
natural to expect the signal to deviate from a simple superposition of linear
QNMs. Such deviations could include new driven frequencies that are
combinations of the first-order QNFs [1389, 1390], or corrections to the first
order modes, in terms of the mode amplitudes, phases, and frequency spectrum
[1391, 1392]. Indeed, a driven second-order mode (which in quasicircular
mergers appears in the $\ell=m=4$ multipole, sourced by the square of the
$\ell=m=2$ fundamental mode) was recently identified in binary BH simulations
[1361, 1377, 1393, 1376, 1394]. Notably, the amplitude of the second-order
mode was comparable to the amplitude from linear order, raising questions
about perturbative convergence.

Detailed second-order calculations remain challenging, and are complicated by
aspects of metric reconstruction and regularization of singularities.
Nevertheless, second-order Kerr perturbations in the ringdown context were
recently calculated numerically [1392]. Additional analytic progress in the
higher-order ringdown has included the development of a bilinear form under
which Kerr QNMs are orthogonal [1395, 1396] (see also [1397]), new methods to
reconstruct the metric in the presence of a first-order source [1305], and the
use of analytic approximations [1390]. These calculations have considerable
overlap with the self-force problem, where significant progress has recently
been made at second order, see Sec. 4.3.3.

####  4.4.4 Challenges

The main challenges in post-merger signal modeling using BHPT concern the
inclusion of higher modes, the instability and quick damping of overtones,
non-QNM ringdown components (e.g., from tails or the prompt response), and
nonlinear effects.

Higher $(\ell,m,n)$ modes are predicted by BHPT to arise in generic binary
mergers. However, there is no consensus on whether it is consistent to include
more than a few overtones ($n\geq 1$) in a purely linear model [1369, 1365,
1380]. The challenge lies in the fact that higher overtones decay rapidly, and
therefore are relevant at early times, when nonlinearities are also expected
to become more significant. Despite providing good fits to the signal, higher
overtones might simply play the role of fitting noise, nonlinearities, or non-
QNM components (such as tails or the prompt response) in the signal. For
asymmetric binaries (with precession, eccentricity, or unequal mass ratios)
angular modes beyond $\ell=m=2$ may also have significant amplitudes (see e.g.
[608, 1354, 1373, 1398]), and their fundamental ($n=0$) modes have lifetimes
similar to the fundamental $\ell=m=2$ mode.

Recent theoretical studies also indicate that the QNM spectrum itself could be
unstable [1399, 1400, 1401, 1402, 1403, 1404, 1405], with overtones
particularly sensitive to nonlinear perturbations, small environmental
effects, and deviations from vacuum general relativity. However the
experimental implications of this instability are not fully understood.

The size of power-law tail contributions due to back-scattering has yet to be
estimated for binary BHs in full GR. Reference [1365] cautioned that, at the
linear level, the tail contribution can be comparable to that of high-overtone
QNMs with $n\simeq 5$. This should motivate a more careful study of the
contribution of back-scattering effects in the post-merger signal.

No study has yet quantified which (if any) LISA-band ringdown sources will
require higher-order perturbation theory. Ideally, starting from a first-order
perturbation, a nonlinear model would predict detailed corrections including
the amplitude of driven modes, shifts in amplitudes and phases for first-order
modes, and any frequency drifts. Some initial progress towards such a model
includes agnostic fits of numerical relativity waveform catalogues [1365],
numerical studies [1392], and developments in BH perturbation theory [1395,
1306].

Some studies have speculated that nonlinearities could be stronger for near-
extreme remnants, whose spectrum contains long-lived modes with commensurate
frequencies. In this limit, nonlinear effects could lead to gravitational
turbulence [1406] and connect with the Aretakis instability of extremal BHs
[1407]. To assess these possibilities, it will be necessary to extend
calculations beyond scalar-field toy models [1406] to full general relativity,
see for example Ref. [1408].

###  4.5 Effective-one-body waveform models

Coordinators: Tanja Hinderer, Geraint Pratten
Contributors: S. Akcay, A. Antonelli, S. Bernuzzi, A. Buonanno, J. Garcia-
Bellido, A. Nagar, L. Pompili

####  4.5.1 General description

The effective-one-body (EOB) approach was originally introduced in Refs. [755,
1100, 1409, 1410, 1244] with the aim of providing GW detectors with
semianalytic waveform models for the entire coalescence of compact-object
binaries (i.e., the inspiral, plunge, merger and ringdown), resumming PN
information around the strong-field test-body limit. Since the breakthrough in
NR in 2005 [514, 515, 516], the EOB framework has incorporated information
from the NR simulations, thus producing highly-accurate waveform models for GW
observations (e.g., see the review articles [1411, 800] and discussion below).
Over the years, the EOB framework has been extended to the scattering problem,
and has incorporated analytical information from other methods, such as the PM