approach and GSF theory, as illustrated in Fig. 9.

The EOB waveform models consist of three main building blocks: 1) the
Hamiltonian, which describes the conservative dynamics, 2) the radiation-
reaction (RR) force, which accounts for the energy and angular momentum losses
due to GW emission, and 3) the inspiral-merger-ringdown waveform modes, built
upon improved PN resummations for the inspiral part, and functional forms
calibrated to NR waveforms for the merger-ringdown signal. We now briefly
review these three key ingredients.

A fundamental pillar of the EOB approach is the map of the real two-body
dynamics into that of an effective test mass or test spin in a deformed
Schwarzschild or Kerr background, with the deformation parameter being the
symmetric mass ratio $\nu=\mu/M$, where $\mu=m_{1}m_{2}/M$ is the reduced mass
and $M=m_{1}+m_{2}$ is the total mass. More specifically, the real or EOB
Hamiltonian, $H_{\mathrm{E}OB}$, is related to the effective Hamiltonian,
$H_{\mathrm{e}ff}$, by [755]

$H_{\mathrm{E}OB}=M\sqrt{1+2\nu\left(\frac{H_{\mathrm{e}ff}}{\mu}-1\right)}.$
(15)

Interestingly, such a relation agrees with calculations in quantum
electrodynamics [1412] aimed at deriving an approximate binding energy for
charged particles with comparable masses in the eikonal approximation; it has
been shown to hold exactly at 1PM order [915], and it has been extensively
used in scattering-theory computations [915, 970, 983, 916]. The above energy-
map (15) achieves a concise resummation of PN information into the Hamiltonian
via a small number of terms. In the center-of-mass frame, the EOB equations of
motion read:

$\frac{d\bm{r}}{dt}=\frac{\partial
H_{\mathrm{E}OB}}{\partial\bm{p}}\,,\quad\quad\frac{d\bm{p}}{dt}=-\frac{\partial
H_{\mathrm{E}OB}}{\partial\bm{r}}+\bm{\mathcal{F}}\,,\hskip
25.0pt\frac{d\bm{S}_{1,2}}{dt}=\frac{\partial
H_{\mathrm{E}OB}}{\partial\bm{S}_{1,2}}\times\bm{S}_{1,2},$ (16)

where $\bm{r}$ and $\bm{p}$ are the canonical variables, notably the relative
position and momentum, respectively, $\bm{\mathcal{F}}$ is the RR force, and
$\bm{S}_{i}$ with $i=1,2$ are the spins of the compact objects. For example,
in the nonspinning limit, where the dynamical variables reduce to
$(r,\phi,p_{r},p_{\phi})$, the effective Hamiltonian in the gauge of Refs.
[755, 1409] reads

$H_{\text{eff}}=\mu\,\sqrt{A_{\nu}(r)\left[\mu^{2}+A_{\nu}(r)\bar{D}_{\nu}(r)\,p_{r}^{2}+\frac{p_{\phi}^{2}}{r^{2}}+Q_{\nu}(r,p_{r})\right]}\,,$
(17)

where the potentials $A_{\nu}(r)$ and $D_{\nu}(r)$ differ from the
Schwarzschild ones due to PN corrections depending on $\nu$. The higher-order
(yet) unknown PN corrections in $A_{\nu}(r)$ might be informed to NR
simulations. The potential $Q_{\nu}(r)$ is a non-geodesic term that needs to
be introduced at 3PN order [1409] to preserve the mapping (15), and reduces to
zero in the test-mass limit.

Precessing-spin EOB waveforms for the inspiral-merger-ringdown were first
built in Ref. [1244] using the nonspinning EOB Hamiltonian [755, 1100, 1409]
augmented with a spinning PN Hamiltonian [767]. The EOB Hamiltonian for
spinning objects was first developed in Ref. [1410], and then in Refs [1413,
1414, 1415, 1416]. These papers follow the structure of the Hamiltonian
advocated in Ref. [1410] that in the large mass-ratio limit reduces to the one
of a (nonspinning) test mass on a Kerr background. Furthermore, another line
of research, which started in Refs. [1417, 1418], built EOB Hamiltonians for
spinning objects such that, in the large mass-ratio limit, they reduce to the
one of a test spin on a Kerr background [1417, 1419]. The two different
spinning Hamiltonians were comprehensively compared in Refs. [1420, 1421].

Figure 9: Theoretical inputs to the EOB framework. The EOB theory draws on a
variety of perturbative results from numerous approaches in different regimes,
as well as from NR simulations, as illustrated here.

In the original EOB model [1100], the radiation-reaction force of Eq. (16) was
given by a suitable Padé resummation of the PN-expanded energy flux at 2.5PN
order, following the seminal work of Ref. [1422]. Subsequently, the (quasi-
circular) radiation-reaction force, that is the flux of angular momentum, has
been expressed as the sum of factorized and resummed multipoles according to
the procedure introduced in Refs. [1246, 1423]. This approach was then first
extended to spinning bodies in Ref. [1158] and then improved in Refs. [1424,
1425, 1426, 1427, 1275, 1428] by means of additional factorizations and
resummation of the orbital and spin parts, and the inclusion of higher-order
PN terms.

The factorization of each waveform multipole proposed in Ref. [1423] reads

$h_{\ell m}^{\text{insp-plunge}}(t)=h_{\ell
m}^{(N,\epsilon)}\,S_{\text{eff}}^{(\epsilon)}\,T_{\ell m}\,e^{i\delta_{\ell
m}}\,(\rho_{\ell m})^{\ell}\,{h}^{\mathrm{N}QC}_{\ell m}\,,$ (18)

where $h_{\ell m}^{(N,\epsilon)}$ is the Newtonian contribution and $\epsilon$
denotes the parity of the mode. The factor $T_{\ell m}$ resums an infinite
number of leading-order logarithms arising from tail effects [1246] (see also
Ref. [826]), the term $e^{i\delta_{\ell m}}$ is a residual phase correction
due to sub-leading order logarithms in hereditary contributions, and the
functions $\rho_{\ell m}$ are the residual amplitude corrections. These
functions were originally obtained as PN-expansions [1423], but in some
versions of the EOB models, it was found useful to further resum them using
Padé approximants for improved strong-field robustness [1426, 1427, 1275].
Another approach consists of calibrating effective high-order PN parameters
entering the $\rho_{\ell m}$’s to improve their accuracy, doing so either
using NR data [1424, 1425, 1428, 1429] or second-order GSF results [1430].

The factor ${h}^{\mathrm{N}QC}_{\ell m}$ is the phenomenological next-to-
quasi-circular (NQC) correction [1431] to the waveform that is informed by NR
simulations, and it is designed to correctly shape the waveform during the
late plunge up to merger, where the motion is not quasi-circular and the
resummed quasi-circular waveform lacks information. The complete EOB waveform
is constructed by attaching the merger-ringdown mode, $h_{\ell
m}^{\mathrm{m}erger-RD}(t)$, to the inspiral-plunge one, $h_{\ell
m}^{\mathrm{i}nsp-plunge}(t)$, at a suitable matching time
$t=t_{\mathrm{m}atch}$, around the peak of the EOB orbital frequency (which
approximately corresponds to merger time), that is [1100]:

$h_{\ell m}(t)=h_{\ell m}^{\mathrm{i}nsp-plunge-
merger}(t)\,\Theta\left(t_{\mathrm{m}atch}^{\ell m}-t\right)+h_{\ell
m}^{\mathrm{r}ingdown}(t)\,\Theta\left(t-t_{\mathrm{m}atch}^{\ell
m}\right)\,,$ (19)

where $\Theta(t)$ is the Heaviside step function. Inspired by results for the
infall of a test mass into a BH [1432] and the close-limit approximation
[1433], several EOB waveform models used a superposition of QNMs for the
dominant mode [1100, 1410, 575, 1434, 1247, 1435, 1126, 1424, 1436, 1251,
1425, 1437]. More recently, an NR-informed fit of the amplitude and phase has
become standard [1438]. This framework is similar in spirit to the rotating
source approximation developed in Ref. [1439], but with significant technical
differences. An important input into the merger-ringdown part of the EOB model
is the mapping between progenitor binary parameters and the mass and spin of
the final remnant. This is needed to determine the frequency of QNMs of the
remnant, and it is obtained from NR results [1440, 599].

There is considerable freedom in modeling and resumming the EOB Hamiltonian,
radiation-reaction force and GW modes, and spin effects. Those different
choices, together with variations in the gauge adopted and deformations of the
potentials in the Kerr spacetime, have led to two main EOB families: SEOBNR
(e.g., see Refs. [754, 1428, 756, 1441]) and TEOBResumS (e.g., see Refs.
[1415, 1442, 1443, 760]), which we discuss in detail in Secs. 4.5.3 and 4.5.4
below. The EOB waveform program advances along two interrelated directions.
The first comprises theoretical advances on the underlying structure and
mapping of analytical results from PN, PM and GSF approaches. The second
involves testing and refining the models using NR information and rendering
them available for use in data analysis. The EOB waveform models are
intrinsically in the time domain (TD), which provides a more direct relation
between source physics and asymptotic radiation than frequency-domain models.
Modeling both the dynamics and radiation is useful for gaining deeper insights
into strong-field nonlinearities and it enables more thorough tests of the
robustness of different aspects of the waveform based on examining the
behavior of various characteristic quantities under changes in the binary’s
parameters, and on comparing different gauge-invariant quantities between EOB,
NR, or perturbative results — for example the binding energy and fluxes [1444,
1445, 1005, 1446, 1447, 917], periastron advance [1114, 1115, 1448],
scattering angles [1449, 720, 986, 987] or spin-precession invariant [1448]
and redshift. At the same time, this also makes EOB models computationally
more expensive than IMR waveform models, which describe only the asymptotic GW
signals in the frequency domain (FD), as discussed in Sec. 4.6. However, the
computational cost is not an insurmountable problem. Much recent work has
significantly improved the computational efficiency of EOB models through the
development of optimized codes [1450, 1451, 1452, 1420, 1453], by building
reduced-order models (ROMs) [1454, 1455, 1456, 1457, 1458, 1459, 1460, 1461,
1462, 1463, 1429], implementing the quasi-adiabatic approximation [1452, 1420,
760, 1464], and making use of machine-learning methods [1465, 1466, 1461,
1463]. Close Hyperbolic Encounters could also be searched for with both
ground- and space-based interferometers using machine learning [1467].

EOB waveform models have been constructed for quasi-circular non-spinning
[1100, 1409, 575, 1434, 1468, 1423, 1469, 1435, 1436, 805, 1470] and spinning
[1410, 1244, 1471, 1413, 1472, 1473, 1418, 1419, 1414, 1415, 1474, 1416, 1421,
1424, 1425, 754, 1428, 1251, 1437, 756, 1475, 1442, 1476, 760, 1429, 1441]
binaries. Furthermore, orbital eccentricity [1477, 752, 714, 759, 757, 758,
1478] and matter effects [1446, 1479, 1480, 1481, 1482, 1483, 1484] (see Sec.
6.5) , as well as information from PM [915, 916, 917, 985, 986, 987] and
conservative GSF information [1485, 1486, 1487, 1009, 1006, 984, 1488, 505,
1489, 1430] have been also incorporated in EOB models. The required effort to
advance the EOB modeling involves: 1) obtaining new analytical information and
mapping it into the EOB framework; 2) resumming the results into full waveform
models with additional flexibility for calibrations to NR; and 3) testing and
comparing models to assess their performance and accuracy. We now describe
each of on the two main EOB families in some more detail.

family | waveform model | $spins$ | ecc./hyper. | NR & Teukolsky calib. region | $CP-framemodes$
---|---|---|---|---|---
$\mskip 5.0mu plus 5.0mu\mskip 5.0mu plus 5.0mu1^{\mathrm{s}t}$ |  ∗EOBNRv1 [1434] |  |  | $q\leq 4$ | $(2,2)$
$\mskip 5.0mu plus 5.0mu\mskip 5.0mu plus 5.0mu2^{\mathrm{n}d}$ |  EOBNRv2 [1469] |  |  | $q\leq 3$ | $(2,2)$
|  SEOBNRv0 [1472] |  |  | $q=1,\chi_{1,2}=\pm 0.4$ | $(2,2)$
|  ∗EOBNRv2HM [1435] |  |  | $q\leq 6$ | $(2,2),(2,1),(3,3),(4,4),(5,5)$
|  ∗EOBNRv2HM_ROM [504] |  |  | $q\leq 6$ | $\\!\\!\\!\\!(2,2),(2,1),(3,3),(4,4),(5,5)\\!\\!\\!\\!$
$\mskip 5.0mu plus 5.0mu\mskip 5.0mu plus 5.0mu3^{\mathrm{r}d}$ |  ∗SEOBNRv1 [1424] | $\checkmark$ |  | $q\leq 6,\chi_{1,2}=\pm 0.4$ | $(2,2)$
|  ∗SEOBNRv1_ROM [1455] | $\checkmark$ |  |  |
|  ∗SEOBNRv2 [1425] | $\checkmark$ |  | NR: $q\leq 8,|\chi_{1,2}|\leq 0.98$ |
Teuk.: $q=1000,|\chi_{1,2}|\leq 0.99$ |  |  |  |  |
|  ∗SEOBNRv2_ROM [1456] | $\checkmark$ |  |  |
|  SEOBNRv2P [1251] | $\checkmark\checkmark$ |  |  | $(2,2),(2,1)$
|  ∗SEOBNRv3P [1425, 1437] | $\checkmark\checkmark$ |  |  |
$\mskip 5.0mu plus 5.0mu\mskip 5.0mu plus 5.0mu4^{\mathrm{t}h}$ |  ∗SEOBNRv4 [754] | $\checkmark$ |  | NR: $q\leq 8,|\chi_{1,2}|\leq 0.995$ |
Teuk.: $q=1000,|\chi_{1,2}|\leq 0.99$ | $(2,2)$ |  |  |  |
|  ∗SEOBNRv4_ROM [754] | $\checkmark$ |  |  | $(2,2)$
|  ∗SEOBNRv4P [756] | $\checkmark\checkmark$ |  |  | $(2,2),(2,1)$
|  ∗SEOBNRv4HM [1428] | $\checkmark$ |  |  | $\\!\\!\\!\\!(2,2),(2,1),(3,3),(4,4),(5,5)\\!\\!\\!\\!$
|  ∗SEOBNRv4HM_ROM [1459] | $\checkmark$ |  |  |
|  ∗SEOBNRv4PHM [756] | $\checkmark\checkmark$ |  |  |
|  ∗SEOBNRv4PHM_surr [1460] | $\checkmark\checkmark$ |  |  |
|  ∗SEOBNRv4EHM [758] | $\checkmark$ | $\checkmark$ |  |
|  ∗SEOBNRv4PHM_4dq2 [1463] | $\checkmark\checkmark$ |  |  | $(2,2),(2,1),(3,3),(4,4)$
$\mskip 5.0mu plus 5.0mu\mskip 5.0mu plus 5.0mu5^{\mathrm{t}h}$ |  ∗SEOBNRv5 [1429] | $\checkmark$ |  | NR: $q\leq 20,|\chi_{1,2}|\leq 0.998$ |
Teuk.: $q=1000,|\chi_{1,2}|\leq 0.99$ | $(2,2)$ |  |  |  |
|  ∗SEOBNRv5_ROM [1429] | $\checkmark$ |  |  |
|  ∗SEOBNRv5HM [1429] | $\checkmark$ |  |  |
(4,3),(4,4),(5,5) |  |  |  |  |
|  ∗SEOBNRv5PHM [1441] | $\checkmark\checkmark$ |  |  |
$\checkmark$ aligned spin, $\checkmark\checkmark$ arbitrary spin orientations. |  |  |  |  |
∗ implemented in LALsuite, and for the latest version (v5) in the open-source Python package pySEOBNR [1490]. |  |  |  |  |

Table 9: Progress in the development of ${\mathtt{S}EOBNR}$ models. The NR-
calibration region refers only to calibration against aligned-spin NR
waveforms. We indicate the coprecessing modes as CP modes. The surrogate
models ${\mathtt{S}EOBNRv4PHM\\_surr}$ and $\mathtt{SEOBNRv4PHM\\_4dq2}$, and
the eccentric/hyperbolic model SEOBNRv4EHM are implemented in LALSuite, but
they have not been reviewed, yet, and they are not publicly available at the
moment. The surrogate models are limited in length and binary’s parameter
space (see Refs. [1460, 1463] for details).

####  4.5.2 Suitable for what sources?

Inspiral-merger-ringdown signals from massive black-hole binaries (MBHBs) will
be within the LISA band (see Sec. 2.1), making EOB models naturally applicable
for these sources. The bulk of MBHB events observable by LISA are expected to
have mass ratios $q\lesssim 10$, but some systems may have mass ratios of up
to several hundreds. Most tests and calibrations of EOB models have been for
$q\lesssim 20$, but the fact that EOB models interpolate between this regime
and the test-body limit with information from BH perturbation theory makes
them structurally well-suited for larger-$q$ systems. Spin precession and
eccentricity effects, both simultaneously relevant for MBHBs and IMRIs, have
to date been separately included in EOB models. These considerations apply not
only for MBHBs but also for describing intermediate mass-ratio binaries
(IMRIs) discussed in Sec. 2.3. Importantly, EOB models have been either
validated [505, 1489] or improved using perturbation-theory and GSF
information [1486, 1487, 1009, 984, 1488, 1430], which is relevant for IMRIs
and also EMRIs. EOB models are also suitable to stellar-origin binaries Sec.
2.5, whose early-inspiral signals will be within the LISA band. For inspirals,
EOB models provide a more accurate though less efficient description than pure
PN-based models, where, e.g., Ref. [1453] shows that using EOBSPA+PA is rather
efficient while searches using PN waveforms would miss a significant portion
of events. Using EOB waveforms will further be important for connecting the
LISA portion of SOBHB signals with the corresponding merger signals measurable
in future ground-based detectors. The EOB approach also allows for the
inclusion of additional physical effects, for instance, from gravity theories
beyond GR and environmental effects as explained in Sec. 6.5.

As discussed in Sec. 4.5.6, while EOB models are naturally suited to the above
sources, significant further advances in the modeling accuracy, complexity of
physical effects, robustness over a wide range of parameter space, and
efficiency will be required for their use in LISA data analysis.

####  4.5.3 The SEOBNR waveform models

The EOBNR family of waveform models 333The generic name SEOBNRvnEPHM indicates
that the version vn of the EOB model is calibrated to NR simulations (NR),
includes spin (S) and precessional (P) effects, eccentricity (E) and higher
modes (HM). has been developed with two main goals: i) make use of the most
accurate analytical information for the two-body dynamics and gravitational
radiation (PN, PM, GSF) and results from NR and Teukolsky-code simulations to
build physical, highly-accurate waveform models of compact-object binaries,
and ii) make them available to GW detectors for searches and inference
studies. The history, since 2007, and the current status of EOBNR models are
illustrated in Table 4.5.1. These models have been implemented into the LIGO
Algorithms Library (LALsuite) [1491], and more recently in the open-source
Python package pySEOBNR [1490]. These codes were then reviewed by the LVK
Collaboration. Here we summarize the main milestones, focusing on the binary
BHs, while leaving details of the modeling improvements (two-body inspiraling
dynamics, transition merger to ringdown, RR effects, GW modes, resummation of
EOB potentials, etc.) to the corresponding publications.

Building on the original work on the EOB framework [755, 1100, 1409, 1244],
the first EOBNR model calibrated to nonspinning NR waveforms was developed in
Ref. [1434], following initial comparisons to NR equal-mass nonspinning
binaries in Ref. [575], in which the importance of including overtones in the
EOB description of the merger-ringdown signal was pointed out. EOBNRv1 was
employed by Initial and Enhanced LIGO, and Virgo for the first searches of
coalescing binary BHs [1492, 1493, 1494, 1495]. The calibration to spinning
binaries with equal masses was considered in Ref. [1472], employing an EOB
Hamiltonian for a test mass in a deformed Kerr spacetime [1413]. The
nonspinning waveforms with higher modes were first modeled in Ref. [1435] with
the improved factorized waveforms [1246, 1423, 1158], thus marking the second
generation of EOBNR models.

The third-generation of EOBNR models encompassed significant advances in
including spin effects in the two-body dynamics and radiation, resumming
perturbative information, and calibration to NR simulations. An EOB
Hamiltonian for a test spin in a deformed Kerr spacetime was derived in Refs.
[1417, 1418, 1419]. It included all PN corrections in the test-body limit, at
linear order in the test-body’s spin. To enforce the presence of a photon
orbit (and peak of the orbital frequency) for aligned-spin binary BHs, it used
a logarithmic (instead of Padé) resummation of the EOB potentials. Those
improvements led to the aligned-spin model SEOBNRv1 [1424], and the first
spin-precessing model, SEOBNRv2P [1251], which adopted the co-precessing frame
description of Refs. [1496, 1497] to efficiently handle precessional effects.
Those models also included information for the merger-ringdown waveforms from
the test-body limit, notably from time-domain Teukolsky waveforms obtained in
Refs. [1126]. With the availability of a larger set of accurate SXS NR
waveforms, and further analytical results from PN theory, the model was
upgraded to SEOBNRv2 [1425], with its precessing version known as SEOBNRv3P
[1437]. These third-generation SEOBNR models were used in the template bank of
the modeled searches (gstLAL and pyCBC), and for parameter-estimation of GW
signals detected during the first observing run (O1) of Advanced LIGO [193,
770].

The fourth generation of models included substantial and more efficient (via
Markov-chain–Monte-Carlo methods) re-calibrations of the aligned-spin baseline
models, such as SEOBNRv4 [754], and phenomenological ansatzes for the merger-
ringdown based on Ref. [1438]. The SEOBNRv4 model was completed by including
higher modes [1428], arbitrary spin orientations [756], and extended to
eccentric and hyperbolic orbits for aligned spins in Ref. [758], using GW
modes recast in factorized form in Ref. [757] (see also Refs. [713, 715, 1498,
752] for other eccentric EOBNR models). These updated models for quasi-
circular orbits were used for the template banks of modeled searches, and
extensively for inference studies during the second (O2) and third (O3) runs
of Advanced LIGO and Virgo [1499, 113, 115, 1500, 210].

Among the highlights of the fifth generation of models, they were calibrated
to NR waveforms with larger mass-ratios and spins using a catalog of 442 SXS
simulations [1429], and they incorporate for the first-time information from
second-order GSF in the nonspinning modes and radiation-reaction force [1430].
The models include several higher modes, notably the $|m|=\ell$ modes for
$\ell=2\mbox{--}5$, the $(2,|1|)$, $(3,|2|)$ and $(4,|3|)$ modes. Differently
from the SEOBNRv3 and SEOBNRv4 Hamiltonians, the SEOBNRv5 reduces in the test-
body limit to the one of a test mass in the Kerr spacetime. The accurate
precessing-spin dynamics of SEOBNRv5PHM [1441] was obtained including partial
precessing-spin information in the EOB Hamiltonian in the co-precessing frame
by orbit averaging the in-plane spin components of the full precessing
Hamiltonian at 4PN order [1501]. Furthermore, the evolution equations for the
spin and angular momentum vectors were computed in a PN-expanded, orbit-
averaged form for quasi-circular orbits, similarly to what was done in Refs.
[1502, 1503, 760], but with important differences due to the distinct gauge
and spin supplementary conditions [1501], and the inclusion of orbit-average
in-plane spin effects. Thanks to a simpler (but more approximated) precessing-
spin inspiraling dynamics, which allows for the use of the post-adiabatic
approximation [1420, 1464], and the new high-performance Python infrastructure
pySEOBNR [1490], the computational efficiency of the SEOBNRv5PHM models has
improved significantly, generally $\sim 8\mbox{--}20$ times faster than the
previous generation SEOBNRv4PHM. This is particular important in view of the
use of these waveform models with next-generation detectors on the ground and
with LISA, which have a much broader bandwidth.

Furthermore, SEOBNR waveform models have also been used to extend the NR
surrogate models (see Sec. 5.1) to lower frequencies, by hybridizing them
[1504, 765, 1505], and to build IMR phenomenological waveform models (see Sec.
4.6). For the latter, SEOBNR models have been employed to construct time-
domain hybrid NR waveforms, and also to calibrate the phenomenological models
in regions of the parameter space where NR waveforms are not available [728,
762, 1506].

Quite importantly for LISA and next-generation detectors on the ground, the
computational efficiency of the time-domain SEOBNR models can also be improved
by building surrogate and ROM versions [1455, 1456, 754, 1460, 1429]. Lastly,
the efficient and flexible pySEOBNR infrastructure [1490], which uses Bayesian
algorithms to calibrate waveforms against NR simulations, will allow to
swiftly include and test new analytical information (PN,PM,GSF) (and their
possible resummations) as soon as they become available. This is crucial to
improve waveform models by at least two orders of magnitude to match the
expected waveform accuracy requirements.

Extensions of the SEOBNR models to extreme mass-ratio inspirals were obtained
in Refs. [1486, 1487, 1009, 1425, 1127], and to non-vacuum binary systems and
gravity theories beyond GR are discussed in Sec. 6.5.

####  4.5.4 The TEOBResumS waveform models

Physical content | GIOTTO | DALI
---|---|---
Analytic information | $A(r)$ | Pseudo $5$PN, resummed | Pseudo $5$PN, resummed
$D(r)$ |  $3$PN, resummed |  $5$PN, resummed
$Q(r,p_{r_{*}})$ |  $3$PN |  $5$PN (local)
$G_{S},G_{S_{*}}$ |  $3.5$PN |  $3.5$PN
$r_{c}$ | NLO | NLO
NR information | $a^{6}_{c}$ | Effective parameter in $A(r)$
$c_{\mathrm{N}3LO}$ | Effective parameter in $G_{S},G_{S_{*}}$
NQCs | Ensure correct transition between plunge and merger
Ringdown | Phenomenological model, quasi-circular
BBH NR Validation region | $q\geq 10$ and test-mass with $|\chi_{1,2}|\geq 0.99$; $10<q<128$ no-spins
Spins | Aligned | $\checkmark$ | $\checkmark$
Precessing | $\checkmark$ | -
Orbital dynamics |  | Circular | Generic (bound & open)
Subdominant CP modes | Inspiral to merger | $(\ell,|m|)\leq 8$
Merger/ringdown | $(\ell,|m|)=(2,1),(2,2),(3,2),(3,3),(4,2),(4,4),(5,5)$
Table 10: Current/default physics incorporated in TEOBResumS. TEOBResumS is
developed open source and publicly available at
https://bitbucket.org/eob_ihes/teobresums/. Symbols are defined in the text.
Historical milestones in the model developments and the associated references
can be found in the Wiki page of the repository. Robustness tests and the
detailed parameter space coverage can be found in the Wiki page and are
continuously updated. TEOBResumS can also be installed via pip install
teobresums. The code is interfaced to state-of-art GW data-analysis pipelines,
including bilby and pycbc. The code uses semantic versioning since the
deployment of the GIOTTO version. Earlier versions of TEOBResumS are
implemented in LAL Simulation and reviewed by LVK.

The TEOBResumS model builds upon early EOB developments at the Institut des
Hautes Etudes Scientifiques (IHES) that define the basics structure of the
model. The name suggests the key features: arbitrary compact binaries with
tidal (T) and generic spins (S) interactions are modeled by making systematic
use of analytical resummations. This name has first appeared in the tidal
model of Ref. [1446] and in the tidal and spinning model of Ref. [1442]. The
latter include the factorized waveform described above [1423, 1247], the
systematic resummation of EOB potentials via Padé approximants [1507, 1413,
1247, 1414, 1436], NQC corrections [1245, 1246, 1471, 1473, 1436], attachment
of the ringdown around the peak of the orbital frequency (i.e., light-ring
crossing) [1246, 1436], the concept of centrifugal radius for incorporating
spin-spin interactions [1415], resummed gyro-gravitomagnetic functions [1413,
1414, 1415] and the factorized NR-informed ringdown waveform [1438]. These
developments made heavy use of results in the test-mass limit obtained by
means of a new approach to black hole perturbation theory [1245, 1246, 1249,
1250, 1128] to understand each physical element entering the structure of the
waveform and striving for physical completeness, simplicity, accuracy and
efficiency. Currently, TEOBResumS-GIOTTO is a unified framework to generate
complete inspiral-merger-postmerger waveforms for any type of quasi-circular
compact binary (BBH, BNS, BHNS)[1508, 1484, 1509]. TEOBResumS-Dalí is the
model’s extension to generic orbits and describes bound orbits with arbitrary
eccentricity [714, 759, 1510, 1511], hyperbolic orbits and scattering [1512,
1513, 720]. The key features of the models are summarized in Table 10 and
described in the following. An extension of the model that relies on GSF-
informed potentials [1010, 1488] so as to generate extreme-mass-ratio-
inspirals (EMRIs) is available [1489, 505], and work on scalar-tensor gravity
(building upon Refs. [1514, 1515]) is also currently in progress [1516, 1517,
1518].

The structure of the spin-aligned effective TEOBResumS Hamiltonian is

$H_{\text{eff}}=H_{\text{eff}}^{\text{orbital}}+H^{\text{SO}}_{\text{eff}},$
(20)

where $H_{\text{eff}}^{\text{orbital}}$ incorporates even-in-spin
contributions (spin-spin couplings), while $H^{\text{SO}}_{\text{eff}}$
incorporates odd-in-spin ones (spin-orbit couplings). The spin-orbit
contribution reads

$H^{\text{SO}}_{\text{eff}}=G_{S}S+G_{S_{*}}S_{*},$ (21)

where $(G_{S},G_{S_{*}})$ are the gyro-gravitomagnetic functions and $S\equiv
S_{1}+S_{2}$, while $S_{*}=m_{2}/m_{1}S_{1}+m_{1}/m_{2}S_{2}$ are useful
combinations of the spins [1519, 1413, 1414]. In the large-mass-ratio limit,
$S$ becomes the spin of the largest black hole, while $S_{*}$ is the spin of
the small body. The orbital contribution $H_{\text{eff}}^{\text{orbital}}$
incorporates the three EOB potentials $(A,D,Q)$. The $A$ function employs the
full 4PN information [805] augmented by an effective 5PN term, parametrized by
the coefficient $a_{6}^{c}$, which is informed by NR simulations. In GIOTTO,
both $(D,Q)$ are kept at 3PN accuracy. In the latest versions of Dalí, the
full 5PN information is being currently explored [1509]. The spin-orbit sector
includes next-to-next-to-leading-order (NNLO) analytical information [1414] in
the (resummed) gyro-gravitomagnetic functions [1415] augmented by an effective
$N^{3}LO$ coefficient $c_{3}$ that is informed by NR simulations. Spin-spin
couplings are incorporated using the concept of centrifugal radius [1415],
$r_{c}$, and starting from NLO accuracy. TEOBResumS can deal with generically
oriented spins, thus incorporating the precession of the orbital plane and the
related modulation of the waveform. Spin-precession in the quasicircular case
is incorporated with an efficient yet accurate hybrid EOB-post-Newtonian
scheme based on the common “twist approach” [1476, 760, 1484].

In GIOTTO, the radiation reaction and multipolar waveform (up to $\ell=m=8$)
are implemented as an upgraded version of the factorized and resummed
procedure introduced in Ref. [1423]; see Refs. [1443, 1427]. In addition, NR
information is used to determine NQC waveform corrections (multipole by
multipole) so as to correctly shape the waveform around merger. The latter is
then attached to a NR-informed, phenomenological description of the ringdown
based on Ref. [1443]. Importantly, NR information is included in such a way to
maintain consistency between waveform and fluxes [1508]. The leading-order
horizon absorption terms are implemented in the model since its early version
[1415]. In Dalí, the radiation reaction on generic orbits is obtained by
incorporating generic Newtonian prefactors in the factorized quasicircular EOB
waveform [714, 759]. This approach, which is extended to include resummed 2PN
terms, has been verified against exact Teukolsky fluxes for highly eccentric
and hyperbolic geodesics and currently provides the best available
representation of the radiation reaction [1520, 1521, 1522]. A new paradigm to
incorporate PN results into EOB, which promises a further boost in
performance, has been recently proposed in Ref. [1478]. The inclusion of
matter effects and application to probing the nature of compact objects is
discussed in Sec. 6.5.

On the algorithmic side, the development of TEOBResumS has introduced two key
analytical acceleration techniques, namely the post-adiabatic (PA)
approximation of the EOB dynamics [1452] and the EOBSPA [1453]. The PA
approximation is an iterative procedure to solve the circularized EOB dynamics
at given radii (or frequency) thus bypassing the need to numerically solve the
system of ODEs. Performances are reported in e.g. Refs [1452, 1475, 1508,
1461] and show that TEOBResumS’s waveform generation times are competitive
with some of the fastest waveforms and machine learning approximants for long
signals (i.e. BNS inspirals), while retaining the original EOB accuracy. The
EOBSPA is a technique to generate multipolar FD EOB waveforms using the
stationary phase approximation. Ref. [1453] illustrates that the EOBSPA is
computationally competitive with current phenomenological and surrogate
models, and can generate (virtually) arbitrarily long and faithful waveforms
up to merger for any BNS. Currently, the EOBSPA is also the only alternative
to PN methods for efficiently generating intermediate-mass binary black hole
inspiral waveforms for LISA [1453].

####  4.5.5 Environmental effects

The EOB approach provides considerable freedom for incorporating yet unknown
effects, such as beyond-GR gravity (see Sec. 6.6), exotic BH-like objects, or
environmental effects. Analyzing the presence of such effects critically
relies on accurate GR waveforms for comparison and benefits from the
techniques developed for GR. While the modeling of astrophysical environments
has not yet been explored directly within the EOB waveforms program, there is
in principle no obstruction to incorporating them. For example, perturbative
and numerical studies of environmental effects are already available (as
discussed in Secs. 4.2, 4.3 and 4.1) and have identified key distinctive
features, which could be mapped into EOB models similarly to what has already
been accomplished for other physical effects. Likewise, existing capabilities
of using EOB-baseline models for parameterized tests of GR could be extended
to include environmental effects in a parameterized, phenomenological way at
the level of waveforms. These two approaches could be useful for different
purposes.

####  4.5.6 Challenges

Employing EOB waveforms for GW signals of LISA sources requires a significant
continued effort to improve the accuracy and include all physical effects.
Importantly, the waveform models must incorporate the simultaneous effects of
high spins at generic orientations, eccentricity, and large mass ratios $q$,
all of which also require many more higher modes. Furthermore, we require
additional flexibility to discover new physics from beyond GR and the Standard
Model of particle physics, as discussed in Sec. 6. In terms of accuracy, the
models must describe high signal-to-noise-ratio events of MBHBs, at more than
an order of magnitude higher than the calibration standards achieved to date,
and they have to be tested for the expected LISA detector response.

The effort to meet these requirements involves substantial further model
developments on the theoretical and computational fronts. On the theoretical
side, advances on the structure of the EOB models to incorporate more
information from perturbative schemes are essential. This effort relies on
inputs from higher-order calculations with more realistic physics from PN, PM,
and GSF to develop robust EOB models for a wide range of the parameter space.
While several proof-of-principle studies have been carried out on including PM
and GSF information, more work is needed to incorporate them in full state-of-
the-art models in a way that can readily be updated as more information
becomes available. Methods for including precessing-spins and moderate
eccentricity effects (where the systems circularize before merger) are in
place but only separately. In principle, the EOB Hamiltonian is fully generic
but other components of the model are not, and require further structural
developments to incorporate larger eccentricities, as well as the simultaneous
effect of precessing spins. Further considerations will also be required on
the sets of dynamical variables used to evolve the inspirals, for example,
action-angle variables versus the canonical Cartesian coordinates, or a
mixture of them. The higher complexity when including both precession and
eccentricity effects is challenging, and requires due care to ensure gauge-
invariant comparisons. As for existing models, insights from the small-mass-
ratio limit will likely prove useful to address this challenge.

On the practical side, significant further work on testing, calibrating, and
optimizing the EOB models against NR results is crucial to attain the accuracy
required for LISA sources discussed in Sec. 3.2. This in turn relies on the
availability of accurate NR waveforms over a wider parameter space, as
detailed in Sec. 4.1. Current EOB models are only calibrated in the aligned-
spin sector. However, this would need to be changed, and the calibration be
extended to precessing spins to achieve the much higher accuracy requirements
of LISA sources.

Furthermore, accurate waveforms for LISA that include eccentricity, precessing
spins, large mass ratio, higher modes, and means to test for new physics have
a highly complex structure characterized by a large number of different
frequencies. Accurately capturing these features significantly slows down the
computations of waveforms and enlarges the dimensionality of the parameter
space. While EOB waveforms are a priori less efficient than closed-form
models, there exist a number of approaches to overcome these shortcomings, as
discussed above in Secs. 4.5.1, 4.5.3 and 4.5.4.

###  4.6 Phenomenological waveform models

Coordinators: Sascha Husa, Maria Haney Contributors: M. Colleoni, M. Hannam,
A. Heffernan, J. Thompson

####  4.6.1 Description

After the binary black hole breakthrough of 2005, a pragmatic approach to
developing waveform models for compact binary coalescence was required. Such
models needed to describe the waveform from inspiral to ringdown, be tuned to
NR, and suit a wide range of GW data analysis applications. Key requirements
for broad applicability were (and still are) computational efficiency and
broad coverage of the parameter space. This has led to several generations of
frequency-domain and more recently also time-domain models, which have been
implemented in the open source LALSuite framework [1491] under the name of
IMRPhenom, and which are constructed in terms of piecewise closed form
expressions for the amplitude and phase of spherical harmonic modes. A careful
choice of the closed-form expressions allows the maximal compression of
information about the waveform into a small number of coefficients that vary
across the parameter space. More recently, fast ODE integration techniques
have also been developed to model the evolution of the component spins in
precessing binaries [1523]. The principal objective of the phenomenological
waveform program has been to deliver waveform models that keep up with
requirements of data analysis applications and are updated to become
increasingly accurate as detectors become more sensitive.

To construct the models, one proceeds in three stages: First, an appropriate
piecewise ansatz is developed for simple functions, e.g. the amplitude and
phase of spherical or spheroidal harmonics. This ansatz is split into regions,
for example corresponding to the inspiral, merger, and ringdown, where
physical insight about the different regions can be exploited. Fig. 10 shows
the three regions used in current models, where in each such region closed-
form expressions are developed to approximate a discrete data set of
calibration waveforms. Future upgrades may increase the number of regions to
increase accuracy. The analytical ansatz attempts to incorporate physical
insight, e.g. regarding perturbative information concerning the inspiral and
ringdown. For the inspiral, an ansatz is typically constructed as a
deformation of a post-Newtonian description. In the ringdown, black hole
perturbation theory can be used to link features in the waveform to the
quasinormal frequencies determined by the final mass and spin of a black hole.
The ansatz is then fitted to each waveform in a calibration data set,
resulting in a set of generalized coefficients for each waveform. This stage
is usually referred to as the “direct fit”. The input data consist of NR data
and perturbative descriptions at low frequency, such as post-Newtonian
expansions or EOB models. Most typically, these types of input data are used
in the form of hybrid waveforms, which are constructed by gluing inspiral
approximants to shorter NR data. Finally, the coefficients are interpolated
across the parameter space, which we will refer to as the “parameter space
fit”. In both stages of fits it is essential to avoid both overfitting and
underfitting.

Waveform model development has been increasingly “data driven”, adapting
phenomenological descriptions to the available data sets, while also using
physical insight and perturbative information. The inspiral descriptions, for
example, extend post-Newtonian expansions, and the black hole ringdown is
formulated in terms of exponentially damped oscillations with complex
frequencies that depend on the dimensionless spin of the remnant. In addition,
approximate maps have been used, e.g. to map the non-precessing waveforms to
precessing ones. Such maps:

  * $\bullet$

Allow to increase the region of the parameter space where the models can be
employed before sufficiently extensive catalogues of NR waveforms are
available.

  * $\bullet$

Give guidance for future more sophisticated models that are calibrated to NR.

Figure 10: Illustration of the regions used in Phenom models exemplified by
the IMRPhenomX model. The plots show the Fourier domain amplitude (left) and
phase derivative (right) of the $\ell=2,m=-2$ spherical harmonic mode.
Vertical lines mark the frequencies of the ISCO for the remnant black hole,
the MECO (minimum energy circular orbit as defined in [1524]) and the ringdown
frequency of the fundamental mode. Figure adapted from [762].

The “phenom” approach sketched above facilitates the development of simple and
robust codes and very rapid evaluation of the waveforms, and is also well
suited to benefit from parallelisation, e.g., through GPUs [503] or similar
hardware. Phenom models have already been used in computationally efficient
parameter estimation studies of MBHB sources for LISA [503]. Initially, only
frequency domain models were developed [1525, 1526, 1527, 1528, 753, 522, 728,
1529, 1530, 1531, 1532, 1533, 1534, 1535, 762, 1536, 1537, 761], since they
are naturally adapted to matched filter techniques carried out in the
frequency domain. More recently, time domain models have also been developed
[1502, 1506, 1523], which can simplify modelling complex phenomena such as
precession or eccentricity.

Four generations of such models have been constructed to date in the frequency
domain, progressing from the $\ell=|m|=2$ mode for quasi-circular non-spinning
binaries to multi-mode precessing waveforms, as will be described in Sec.
4.6.3. The flexibility of the phenomenological approach is also well suited to
parameterize unknown information, e.g., beyond-GR effects (see Sec. 6.6).

####  4.6.2 Suitable for what sources?

Massive black hole binaries have already been modelled with the phenom
approach, i.e., comparable-mass black-hole binaries, described in Sec. 2.1.
The bulk of such binaries detected by LISA is expected to have mass ratios
$q=m_{1}/m_{2}$ up to $\sim 10$ (see Sec. 2.1.2), where phenom models have
already been calibrated to NR. The main challenge to model such binaries for
LISA is to significantly increase the accuracy, corresponding to the much
higher expected values for the signal-to-noise-ratio (SNR) with LISA. However,
a tail down to mass ratios of a few times $10^{-3}$ also has to be expected
(Sec. 2.1.2). The latest generation of phenom models has indeed already been
calibrated to perturbative numerical waveforms that correspond to solutions of
the Teukolsky equation [762, 1536] at mass ratios of $\approx 10^{-3}$.
Significant progress will still be required regarding both the input
calibration waveforms and the models to increase the accuracy, in particular
for larger spins and the inspiral. Developments and challenges toward
computing large mass ratio NR waveforms are discussed in Sec. 4.1, including
synergies between traditional numerical methods, analytical and self-force
methods. In principle, the phenom approach may be well suited to some EMRI
systems as well. Extension to EMRIs seems most feasible for non-precessing
quasi-circular systems, but may not be practical for precessing or eccentric
systems, where the signal can be extremely complex.

Stellar-origin binaries (Sec. 2.5) and Galactic binaries (Sec. 2.4) are types
of sources for which the phenomenological approach is likely to provide a well
suited framework. It can be adapted to produce computationally efficient
models that can be tuned to data analysis requirements. However, to our
knowledge, no work has yet been carried out in this direction. Finally, one of
the strengths of the phenom approach is precisely that it allows for the easy
inclusion of approximate phenomenological models of physical effects. In this
sense, parameterised modifications to GR or models of environmental effects
would in principle be straightforward to include.

####  4.6.3 Status

To date four generations of FD models have been developed, as well as a first
generation of time domain models, which correspond to the most recent FD
models in terms of parameter space coverage and accuracy. These models have
been implemented in the LIGO Algorithms library (LALsuite) [1491] and reviewed
by the LIGO-Virgo-KAGRA collaboration. To date, the models have only been
calibrated to non-precessing quasi-circular numerical waveforms. Precession
has been implemented via the “twisting” approximation [753, 1529, 1530, 1531,
1532, 761, 1502], which is based on the fact that in the inspiral the
precession timescale is much smaller than the orbital timescale. Consequently,
the precessing motion contributes little to the GW luminosity, and the
inspiral rate is dominated by the non-precessing dynamics. One can thus
approximately map non-precessing to precessing waveforms, and in a further
step one can use fast PN approximants to describe the Euler angles that rotate
the orbital plane and waveform. An effective analytical single-spin
description for the Euler angles including next-to-next-to-leading order
(NNLO) orbit-averaged PN effects has been developed in [1538], and a double
spin approximation based on a multiple-scale analysis (MSA) has been developed
in [497]. For the third generation of these models, an approximate map has
also been used to construct subdominant harmonics from the $l=|m|=2$ mode, but
this approximation has been dropped in the fourth generation, calibrating
instead to NR. A key shortcoming of implementations of the twisting-up
approach in the frequency domain [753, 1529, 1530, 1531, 1532, 761] has been
the use of the SPA, which is not appropriate for the late inspiral, merger and
ringdown. Possible strategies for going beyond the SPA have been discussed in
[1539], but will require further study and testing. Further shortcomings of
the twisting-up approach have recently been discussed in [1540].

Before discussing the current generation of models we briefly summarize the
historical development of the first three generations:

  1. (i)

The first generation consists of the non-spinning PhenomA model for the
$l=|m|=2$ modes of non-precessing binaries [1525, 1526].

  2. (ii)

Second generation models PhenomB [1527] and PhenomC [1528] included non-
precessing spins, where the two spin degrees of freedom are described by a
single effective spin. Also, the precessing PhenomP [753] model was
constructed from PhenomC using an approximate “twisting-up” map between
precessing and non-precessing waveforms.

  3. (iii)

The PhenomD [522, 728] model for the $l=|m|=2$-mode significantly improved the
fidelity of the phenomenological ansatz and the quality of the NR calibration.
Several models have been derived from IMRPhenomD:

     * $\bullet$

PhenomPv2 [1529] updates the PhenomP model for precession.

     * $\bullet$

PhenomHM [1530] adds subdominant harmonics based on an approximate map from
the $l=|m|=2$ mode to subdominant harmonics.

     * $\bullet$

PhenomPv3 [1531] upgrades the NNLO effective single-spin description of
precession that had been used for the first two versions of PhenomP to a
double-spin description based on the MSA approximation [497].

     * $\bullet$

PhenomPv3HM [1532] includes subdominant harmonics in the twisting-up
construction. Furthermore, matter effects have been included in terms of tidal
phase corrections [1533, 1534] and amplitude corrections of NS-BH systems
[1535].

The PhenomD model has been implemented for use in the LISA data challenge
infrastructure [1541], and has been used in the “Radler” edition of the LISA
data challenge. The PhenomD and PhenomPv2 models have been used for parameter
estimation for observed GW events since the first detection of gravitational
waves, the GW150914 event [1542]. Its derivatives including subdominant
harmonics have also been used to analyse events during the third observing
run. A detailed study of waveform systematics [593] found these models to be
accurate enough for the GW150914 event, but also highlighted the need for
further improvements. This is consistent with a more recent study concerning
the need for further improvements for upcoming upgrades to ground-based GW
detectors [1].

The PhenomX [762, 1536, 1537, 761] family provides a thorough upgrade of third
generation phenom models. PhenomXAS [762] replaces PhenomD for the dominant
mode, increasing the accuracy by roughly two orders of magnitude in terms of
mismatch:

  * $\bullet$

The number of NR waveforms used for calibration is increased from 20 to $\sim
400$.

  * $\bullet$

Teukolsky waveforms are included in the data set in order to extend the model
to EMRI waveforms.

  * $\bullet$

Final spin and mass include EMRI information derived from circular geodesics.

  * $\bullet$

Both spin degrees of freedom are calibrated to numerical data.

  * $\bullet$

Heuristic parameter space fits have been replaced by the systematic
hierarchical fitting approach described in [1440]. An ansatz is selected among
classes of polynomial and rational functions by minimizing not only the RMS
errors, but also information criteria that approximate a full Bayesian
approach, which avoids overfitting and underfitting. This method has, however,
not yet been applied to dimensions larger than 3 (i.e. precession or spinning
eccentric binaries).

PhenomXHM [1536] extends PhenomXAS to the sub-dominant modes, and PhenomXPHM
[761] adds precession via the twisting-up procedure, and allows to switch
between the NNLO and MSA descriptions for the Euler angles. A first
calibration of Phenom models to precessing NR simulations has recently been
presented [763, 1355], paving the way for further increasing the accuracy of
precessing and eventually generic models in the future. Extensions to BNS and
BH-NS systems following the approaches of third generation models are in
progress.

The fourth generation PhenomX family demonstrates that a significant increase
in the accuracy of phenom models can be paired with a significant decrease in
computational cost. Consequently, PhenomX provides the fastest Inspiral-
Merger-Ringdown (IMR) waveforms currently available for data analysis without
GPU acceleration [1537, 761]. Parameter estimation application with GPU
acceleration for PhenomHM has been presented in [503]. Part of the improvement
in efficiency is due to the “multibanding method”, which is based on [1543]:
analytical error estimates determine the grid spacing for interpolation from a
coarse grid, and a standard iterative scheme is used to rapidly evaluate the
complex exponentials required to compute the waveform from the amplitude,
phase, and the Euler angles used in “twisting up”.

Most recently, the time-domain PhenomT model family [1502, 1506, 1523] has
been constructed to mirror the features of PhenomX in the time domain. A key
motivation for the development of time domain models is that they do not
require an analytical approximation to the Fourier transform in order to
obtain explicit expressions for a “twisted-up” precessing waveform. This also
simplifies the incorporation of analytical results, e.g., concerning the
precessing ringdown frequencies; similar simplifications are hoped for when
modelling eccentricity. For future developments, it is foreseen that both
frequency and time domain models will be upgraded in parallel, and that each
of the two “branches” will benefit from progress with the other.

The historical development and current status of phenom models is sketched in
Fig. 11.

Table 11:  Progress in the development of phenomenological waveform models in
frequency domain (FD) and time domain (TD).

####  4.6.4 Environmental effects

The phenom approach provides significant freedom to incorporate additional
features, such as poorly known subdominant GR effects, beyond-GR effects (see
Sec. 6.6), and environmental effects. We note that analysing the presence of
such effects also requires accurate GR waveforms for comparison and will
benefit from the techniques developed for GR. Work toward tests of GR or the
presence of environmental effects may call for improving the accuracy or other
features of GR models. While the modelling of environmental effects has not
yet been explored in the phenomenological waveform program, the
phenomenological models are primed for their addition in two ways:

  1. (i)

The modular structure of phenomenological waveforms should facilitate the
incorporation of known environmental effects, e.g. by augmenting PN
information about the inspiral with information about environmental effects.

  2. (ii)

Environmental effects for which a quantitative model is not yet available
could be incorporated into the phenom ansatz in a parameterised way, similar
to existing GW tests for theory-agnostic deviations from GR.

In both cases it will be useful to have both frequency and time domain
phenomenological waveform models available, allowing one to choose the natural
domain for a given effect.

####  4.6.5 Challenges

Waveform modelling for comparable mass binaries in general faces four main
challenges:

  1. (i)

The availability of a sufficient number of high-accuracy numerical relativity
waveforms throughout the parameter space, including precession and
eccentricity. (This is not specific for phenom models.)

  2. (ii)

The lack of accurate analytical descriptions, in particular for precession and
eccentricity. Here, specific challenges arise for the phenom approach, which
uses “deformed” PN expressions as the basis for constructing a computationally
efficient inspiral ansatz. This will benefit in particular from further
analytical developments in precession and eccentricity.

  3. (iii)

The development of modelling techniques that can produce efficient models from
high-dimensional data sets. For the phenom approach, it may turn out to be
difficult to develop closed form expressions that accurately represent the
full morphology of precessing and eccentric binaries, and to accurately
interpolate the full parameter space without compromising computational speed.

  4. (iv)

Development of an overall data analysis strategy and concrete code framework.
The flexibility and modularity of the phenom approach could be exploited to
develop variants of models with different tradeoffs between accuracy and
speed, or accuracy and broad coverage of parameter space, which can then be
utilised to optimise computational efficiency.

Significant coupling between these challenges is foreseen, e.g. the number and
quality of NR waveforms required will also depend on the progress with
analytical results. In addition, there will be more sophisticated data
analysis approaches, which could exploit the advantages and disadvantages of
different models, or tunable parameters, which could allow one to trade
accuracy for speed. These could eventually give guidance for the development
of waveform models that are designed as an integral part of the data analysis
strategy. These challenges are not specific to LISA, but apply in general to
further model improvements. LISA, however, poses particular challenges due to
the extreme accuracies required for the loudest events.

The detector response itself is significantly more complicated than for
ground-based detectors (see e.g. [1539]), and the problem of developing fast
and accurate waveform models, as well as strategies of how to best employ them
in data analysis applications, resolving tradeoffs of accuracy versus
computational cost, can not be completely decoupled from considerations
regarding the detector response. In the context of GW data analysis for the
current generation of ground based detectors, a software framework for
waveforms has been established [1491], along with best practices for waveform
review. These practices will need to be extended to the context of LISA
waveforms, the more complicated detector response, and the integration into
the LISA data analysis infrastructure. In order to understand these questions,
in particular the systematic errors and their consequences for both searches
and (especially) parameter estimation, it will be necessary that efficient and
well documented implementations of the LISA response functions are publicly
available, integrated into software packages for matched filter calculations
and Bayesian inference.

The non-precessing quasi-circular sector will provide important initial
guidance on how to connect with EMRI descriptions (or to which degree this is
possible). This will provide an arena for toy-model explorations of highly
accurate models in limited regions of parameter space, determining the limits
of accuracy and computational efficiency that can be achieved, as well as
studying issues of accuracy and efficiency related to the LISA detector
response. It is therefore important to continue pushing to higher accuracies
with these models.

The development of accurate phenomenological models of precessing systems is
expected in both the time and frequency domains. Time-domain models will build
on [1502], which is foreseen to serve as a testing ground to calibrate to NR
and describe effects such as the mode asymmetry associated with large recoils
[1544, 1540, 1545]. Further progress in the precessing inspiral will crucially
depend on progress with approximations, such as those based on the MSA [497]
or dynamical renormalisation group (DRG) [1546, 1547]. In order to accurately
describe the late inspiral, merger and ringdown of precessing systems in the
frequency domain, further tests and implementations of the strategies
discussed in [1539] will be required, such as analytical treatments that go
beyond the SPA.

Eccentric waveform modes are again likely to require cross-pollination between
frequency and time domain models (motivating in part the development of
PhenomT in parallel to PhenomX), and significant advances in the development
of analytical descriptions for the inspiral phase of spinning and eccentric
binaries. This is particularly true where binary systems show both precession
and eccentricity. It will be important to understand whether simple approaches
twisting up eccentric waveforms will be sufficient for moderate SNR. Further
investigations will determine if these prescriptions are adequate to determine
binary parameters with sufficient accuracy, such that only a moderate number
of NR waveforms is required to develop an accurate local model for high SNR
events.

A key challenge is to develop models that are accurate across a large region
of parameter space, describing well the changes of waveform morphology, e.g.,
as the mass ratio increases, or as spin alignment migrates between aligned and
anti-aligned with the orbital angular momentum. For current ground-based
detectors, parameter estimation posteriors are typically rather broad, so
refined models in a smaller region of the parameter space would have very
limited use. For LISA, however, it is foreseen that a hierarchy of models will
be developed: those that are sufficiently accurate to identify signals across
large portions of parameter space, and others that further increase the
accuracy for smaller regions of parameter space, which would be used for the
most accurate parameter estimation of the loudest signals. Models with
parameters allowing one to trade accuracy for computational efficiency would
be appropriate for hierarchical approaches to data analysis. This approach
could also be extended further by mixing descriptions in terms of closed-form
expressions with equation-solving, e.g. to improve the accuracy of describing
eccentricity or precession. An important open question is how to incorporate
results from the self-force program (see Sec. 4.3) into phenomenological
models, and how to use ideas developed for the phenomenological waveforms
program to accelerate self-force waveforms.

##  5 Waveform generation acceleration

Coordinator: Alvin Chua, Michael Katz Contributors: S. Field

The stringent efficiency requirements described in Sec. 3.3 naturally
necessitate the acceleration of waveform models, in order to enable the LISA
science analysis. The two main strategies to achieve this are the development
of improved computational techniques, and hardware acceleration.

###  5.1 Computational techniques

To maximize the science returns from highly relativistic LISA sources with
exacting modeling requirements, we require waveform-acceleration techniques
that do not compromise waveform accuracy. The general solution to this problem
is to construct interpolants or fits for high-fidelity waveform data, where
the data comes from an accurate underlying model that is too slow to use
directly in data analysis studies (e.g., as it involves solving PDEs). This
data can describe either the full time-/frequency-series representation of the
waveform, or specific waveform components that are more computationally
expensive. While many waveform acceleration techniques have appeared in the
literature, most rely on the three steps summarized below.

Figure 11:  Waveform acceleration can be achieved by building a surrogate
model through a sequence of three steps. (Step 1) The red solid lines show the
waveform training data evaluated at specific parameter values. (Step 2)
Reduced-order modeling is used to compress the waveform data through a
dimensional-reduction step, and the blue dots represent the data that needs to
be modeled after compression. (Step 3) The blue lines indicate fits for the
compressed waveform data across the parameter space. The yellow dot shows a
generic parameter value, say $\pmb{\lambda}_{0}$, that is not in the training
set. To compute the waveform for this value, each fit is evaluated at
$\pmb{\lambda}_{0}$ (the yellow diamonds), and then the reduced-order
representation is used to reconstruct the decompressed waveform (the dotted
black line).

####  5.1.1 Step 1: Data representation

A gravitational waveform’s parametric dependence is generally too complicated
to model directly with an acceptable degree of accuracy. Instead, the waveform
is typically decomposed into waveform data pieces that are simpler, more
slowly varying functions of parameters and time. If we were to take a machine-
learning viewpoint of the problem, data pieces are features of the data, and
discovering what problem-specific features to extract is a feature engineering
task. How should we define our features?

For non-precessing BBHs in quasi-circular orbits, an amplitude and phase
decomposition for each $(\ell,m)$ angular-harmonic mode is an obvious choice
for our data pieces. For precessing systems, the problem becomes significantly
more complicated as the modes have a rich signal morphology. Here, the general
approach has been to transform the modes to a co-rotating frame where the
system has significantly reduced dynamical features [1548, 589, 590]. It is
often advantageous to further decompose the data into symmetric and asymmetric
mode combinations [1548, 589, 590]. The price of this simplification is that
the frame dynamics must also be modeled, so that the co-rotating modes can be
transformed back into the inertial frame of the detector. For eccentric BBH
systems, which will be important for LISA science, the optimal data
representation is currently unknown.

The highly disparate masses and large separation of timescales in EMRI systems
allow their inspirals to be modeled within the two-timescale framework [1092].
This allows the trajectory of the inspiral through the orbital configuration
phase space to be computed in milliseconds [1110, 1119]. The waveform is then
computed by sampling the metric amplitudes along this inspiral typically by
evaluating a precomputed interpolant of the amplitudes. Most often these
amplitudes are computed in the frequency domain so the waveform is decomposed
into harmonics of not just angular dependence, but also the fundamental
frequencies in the system [1156, 1142]. Generically precessing and eccentric
EMRI waveforms are extremely complex, but the above description effectively
breaks the waveform down into slowly evolving sinusoids. The only drawback is
that there are far more modes to compute and sum for a typical waveform
($\gtrsim 10^{5}$), though in practice many of these do not carry much power
[1110].

####  5.1.2 Step 2: Reduced-order modeling

The waveform quantities of interest often reside in a representation space of
extremely high dimensionality $D$, and must generally be interpolated or
fitted over a parameter space of modest dimensionality $d$ as well. For
example, the full waveform itself for strong-field LISA sources typically has
$D\gtrsim 10^{4}$ and $d\sim 10$. These computational difficulties can often
be mitigated through the use of dimensional-reduction tools [1549, 1550, 1551,
1552, 1553, 1554, 1555, 1556, 1557, 1558, 1559, 1560], the collection of which
is broadly referred to as ROMs.

ROMs aim to remove redundancies in large data sets, making them more amenable
to analysis, and identifying relevant features for further approximation. It
typically begins with a greedy algorithm or singular value decomposition that
identifies the set of $n$ “most important” values in parameter space, whose
associated waveform quantities span a reduced-basis space within which the
model is well approximated [1559, 1560, 1454, 1455]. This serves as both a
dimensional-reduction step as well as a parametric-sampling step.

While invaluable for reducing computational cost on the representation end
($n\ll D$), the compression provided by ROMs does not actually lessen the
intrinsic complexity of the model. Indeed, the successful application of ROMs
relies on sufficient coverage of the parameter space with training data in the
first place—or, more realistically, on truncation of the parameter domain to
compensate. ROM themselves also do not address the issue of interpolating or
fitting the reduced representation at high accuracy, which can be challenging
even in parameter-spaces with $d\gtrsim 3$.

####  5.1.3 Step 3: Interpolation and fitting

The last step is to interpolate or fit each compressed data piece (viewed as
scalar- or vector-field data) over $d>3$ dimensions, and to do so both
accurately ($\lesssim 10^{-6}$ error) and efficiently ($\lesssim 1$ seconds).
In the context of sources and models for ground-based observing, a variety of
choices have appeared, including polynomials [1454], splines [1455], Gaussian
process regression [1561], deep neural networks [1562, 1462, 1463, 1563], and
forward stepwise greedy regression using a custom basis [589]. A recent review
article provides informative comparisons between some methods [1564]. The
accuracy and efficiency goals of models for ground-based observing are more
modest, however, and it is unclear if these methods will continue to work for
LISA MBHs [1, 2], or for EMRIs in full generality [1160]. As such, other
interpolation and regression methods should continue to be explored.

####  5.1.4 Status of techniques

In the GW literature, the term surrogate modeling is sometimes used as a
generalization of ROMs (for the full waveform) to include additional steps
beyond dimensional reduction. ROM surrogates have been built for many
different waveform families, regions of parameter space, and signal durations.
Early models focused on closed-form waveform families largely as a proof-of-
principle exercise [1560, 1559, 1565, 1558, 1557, 1556]. Since then, ROM
techniques have been further developed and refined for various LVK (hence
comparable-mass) applications: EOB [1454, 1455, 1456, 1459, 1460], NR [764,
1548, 589, 590, 1504], multiple subdominant modes [1454, 1459, 764, 1548, 589,
590, 1504], frequency-domain waveforms [1455, 1456, 1459], time-domain
waveforms [1454, 764, 1548, 589, 590, 1504], generically precessing systems
[764, 1548, 589, 590], neutron star inspirals with tidal effects [1457, 1458,
1566], and waveforms with $\sim 10^{4}$ cycles [1455, 1456, 1459, 1504].

When applied to computationally expensive models such as EOB or NR, ROM
surrogates can accelerate the generation of a single waveform by factors of
$10^{2}$ (EOB models; described by ODEs) to $10^{8}$ (NR models; described by
partial differential equations), while being nearly indistinguishable to the
underlying model. EOB surrogates are extensively used as part of LVK
parameter-estimation efforts as well as template-bank detection [770, 1542].
NR surrogates have been used in numerous targeted follow-up studies on
specific BBH events [1500, 114, 1567, 1568, 1569]. As such, ROM techniques
have been essential to the widespread use of both EOB and NR waveforms for
realistic data analysis efforts with LVK data.

The present surrogate framework is ill-suited to EMRI waveform modeling, where
the increased duration and mode complexity combine to increase the information
content of the model by at least 20 orders of magnitude [487, 488, 215].
Surrogates might be useful in special cases of the IMRI regime, where these
difficulties are lessened. There has been recent work on the first models
[1117, 766] for quasi-circular, non-spinning IMRI systems, based on solutions
to the time-domain Teukolsky equation rescaled to fit NR results. The latest
surrogate waveform model in this series covers the last $30,000$M of the
inspiral-merger-ringdown signal, with mass ratios ranging from 2.5 to
$10^{4}$. To put this in context, for a million-Solar-mass system, this is
equivalent to approximately 2 days of signal, which is inadequate for the
expected observable duration of months to years for an EMRI of that mass.

In the angular- and frequency-based decomposition of an EMRI model, the
problem separates naturally into two main parts: that of generating the
inspiral trajectory (from which the waveform phasing is straightforwardly
derived), and that of generating the waveform amplitude for each of the $\sim
10^{5}$ constituent modes. Fast frameworks for EMRI trajectory generation
already exit [747, 1092]. They present no computational difficulties at
adiabatic order (for search); however, their extension to post-adiabatic order
(for inference) will require interpolating the first-order fluxes to a
relative precision of the mass ratio ($\sim\epsilon$) [1217]. Trajectory-level
interpolation with ROM is unlikely to be of practical use, due to issues of
accuracy and parameter-space extensiveness. The accuracy requirements on the
waveform amplitudes are far less stringent, though, and it is here that ROM
compression and neural-network interpolation has been applied to good effect
(in the mode direction, rather than the time/frequency direction) [1160].

####  5.1.5 Future challenges

LISA data analysis has three distinctive features that will require new ideas
in waveform acceleration. First, the waveforms are significantly longer in
duration, which will stress both computational resources and current ROM
techniques. Second, as compared to LVK waveform modeling, LISA accuracy
requirements are significantly more demanding. The fitting schemes will thus
need to deliver higher accuracies while also contending with increased
dimensionality. Finally, the orbital motion of LISA introduces a time-
dependency in the response of the detector, which also needs to be modeled
efficiently for data analysis applications; see, e.g., [1539].

In terms of source-specific challenges, eccentricity is expected to play a
more important role in the modeling and interpretation of LISA MBHs. This will
necessarily increase the problem from 7 to 9 parametric dimensions. While this
should not pose significant difficulty for the dimensional-reduction step
[1560, 1570], the parameter-fitting step will be vastly more difficult, and it
is unlikely that any of the existing techniques will continue to perform well.
For EMRIs, the main difficulty is ensuring that the waveform-amplitude fit
retains its accuracy when extended from eccentric Schwarzschild orbits at
present [1160] to generic Kerr orbits with a much larger parameter and
representation space.

###  5.2 Hardware accelerators / configurations

In addition to computational techniques, generating waveforms with a variety
of hardware accelerators can vastly improve analysis time. Hardware
acceleration, in a general sense, means hardware designed for a specific type
of optimization when compared to general-purpose CPUs. Some examples of
commonly used hardware accelerators include GPUs, Field-Programmable Gate
Arrays (FPGAs), Tensor Processing Units (TPUs), and Artificial Intelligence
(AI)-accelerators.

GPUs are designed specifically for parallel computations. These devices
leverage a large number of compute cores ($\sim 1000s$) and specialized memory
structure to make a large number of independent calculations simultaneously.
From an algorithmic standpoint, many of the methods discussed in Sec. 5.1 can
be implemented on GPUs to achieve considerable gains in efficiency. These
include, but are not limited to, artificial neural networks (matrix
multiplication), interpolation, parallel sorting, and basis transformation.
Additionally, some accurate waveform calculations can be implemented in
parallel, allowing for all or parts of waveforms to be implemented entirely
for GPUs. Programming for GPUs is generally performed in C/C++ or, as of more
recently, Python.

While GPUs take advantage of their parallel architecture, they are still
programmed using specific software for a fixed hardware unit. FPGAs, on the
other hand, represent “programmable hardware.” These devices are programmed at
the device level to implement customized hardware flow designed for specific
tasks. Their most common usage is in neural networks and signal processing.
Programming on FPGAs is done in Verilog or VHDL, although some wrappers in
C/C++ are starting to become available.

####  5.2.1 Status

In recent years, GPUs have seen increased interest from the community for use
in waveform creation, as well as GW-related analysis in general.
Unfortunately, FPGAs are yet to be utilized in the GW community. GPUs have
been used for precessing waveforms, as well as population analysis related to
ground-based observations from the LVK observing network [1571]. In this work,
the waveform and the detector response were calculated directly on the GPU
hardware. To limit the cost of Bayesian inference in ground-based observing,
GPUs were recently used to greatly accelerate the creation of surrogate
waveform models using techniques such as artificial neural networks [1462].
The creation of MBHB waveforms has also been implemented for LISA, including
both the source-frame waveforms and the LISA response function for use in LISA
parameter estimation studies [503] ($\sim$500x acceleration). The generation
scheme for this work calculated, in parallel, accurate waveforms and LISA
response information on sparse grids. Then, to scale to the full data stream,
cubic spline interpolation was leveraged. This technique is generally well-
suited for parallelization.

For EMRIs, GPUs have been used to improve the efficiency of black hole
perturbation theory calculations within the GSF community [1572, 1573, 1574].
In this work, the many harmonics of EMRI metric perturbation are computed in
parallel as the orbit evolves. In terms of data-analysis related waveform
creation for EMRIs, the Augmented Analytic Kludge (AAK) [488] waveform from
the EMRI Kludge Suite was accelerated with GPUs (gpuAAK) for use in the LISA
Data Challenges ($\sim$1000x). Additionally, the first fully relativistic fast
and accurate EMRI waveform that goes beyond the kludge approximations was
implemented on GPUs resulting in a considerable increase in efficiency
($\sim$2000x) [1160, 1110]. Both the gpuAAK and the fast and accurate EMRI
waveforms generate sparse orbital and phase trajectories in serial on a CPU.
These sparse trajectories are then scaled with spline interpolation to the
full waveform sampling rate. This scaling is performed for each separate
harmonic mode so that they can be properly and accurately combined into the
final waveform. It is actually this summation over $\sim 10^{6}-10^{7}$ time
points, each with a combination of $\sim$1000s of harmonic modes, that most
requires the GPU architecture for optimal results.

####  5.2.2 Challenges

A central challenge when working with GPUs, FPGAs, or any other specialized
hardware is the trade-off between efficiency and streamlined code creation and
maintenance. For optimal efficiency, it is generally necessary to implement
code specific to each type of hardware. This could mean implementation in an
entirely different language with, for example, an FPGA, or an adjustment for
different hardware or software requirements for each application. An example
of this is the use of shared memory in GPUs. Shared memory is a small amount
of memory available on the GPU chip, which is much faster to access in
comparison to the larger global memory located off of the chip. Leveraging
shared memory is essential to achieving optimal efficiency for many
applications. There is no use of this form of shared memory on CPUs.
Therefore, it may be desirable to implement a fast GPU version and a mirror
CPU version that effectively uses the same exact code as the GPU version with
slight adjustments. On the other hand, a CPU version optimized for a CPU may
be needed, creating two separate codes. In general, it is helpful to minimize
the amount of code duplication as much as possible, but that may be governed
by efficiency requirements.

Many applications require some parallelization over CPUs in addition to
specific hardware considerations. Understanding how to best leverage all
resources is key to achieving maximal efficiency. With that said, maximal
efficiency should be achieved while maintaining a constant user interface.
This allows for the application of a waveform within a desired program while
using a variety of hardware.

Another general concern about the use of specialized hardware is just that: it
is specialized. Purchasing specialized hardware in place of CPUs will penalize
the breadth at which the hardware can be used. This has generally hurt the
availability of accelerators at large. Similarly, accelerators come at a
variety of prices. Academic-grade GPUs (indicating they have strong
performance in double precision) are generally expensive. However, the newest
generations of accelerators combined with modern algorithms within the GW
field have proven that efficiency is not only highly desirable, but cost
efficient when compared to the equivalent number of CPUs (up to $\sim$ 1000s).
Fortunately, researchers now often have access to clusters of CPUs, GPUs,
TPUs, and AI-accelerators through local, national, and international scale
computing centres. Ultimately, the combination of a variety of hardware will
be necessary to accomplish the scientific goals of the LISA mission.

##  6 Modelling for beyond GR, beyond Standard Model, and cosmic string
sources

Coordinators: Richard Brito and Daniele Vernieri Contributors: J.
Aurrekoetxea, L. Bernard, G. Bozzola, A. Cárdenas-Avendaño, D. Chernoff, K.
Clough, T. Helfer, T. Hinderer, E. Lim, G. Lukes-Gerakopoulos, E. Maggio, A.
Maselli, D. Nichols, J. Novák, M. Okounkova, P. Pani, G. Pappas, V.
Paschalidis, M. Ruiz, A. Toubiana, A. Tsokaros, J. Wachter, B. Wardell, H.
Witek

As explained in Sec. 2.7, the observation of gravitational waves with LISA has
an enormous potential to probe physics beyond GR and the Standard Model of
particle physics, in regimes that have so far not been explored [99, 323]. In
order to fulfill this potential, we need to construct accurate waveforms that
take into account signatures of new physics. As we will discuss in this
Section, this is in fact a very challenging task that for most cases is still
at its early stages of development.

In this section we provide a brief summary of how beyond GR theories and the
effect of beyond Standard Model physics can be modelled using the techniques
presented in Sec. 4. The discussion will be mostly concerned with modelling
the _generation_ of gravitational waves by LISA sources in the presence of new
physics. Beyond-GR effects can also affect the way gravitational waves
_propagate_ on cosmological scales, see e.g. [1575, 1576, 1577, 1578, 1579,
1580, 1581, 1582, 275]. In general, separation of scales allows to add
modifications to the GW propagation on top of any waveform model for which the
generation problem is understood and therefore we do not discuss this
possibility here. More details on this problem can be found in Refs. [1578,
275]. We will also not discuss in detail the potential implications that LISA
observations will have for fundamental physics, since a more detailed review
on this issue can be found in LISA’s Fundamental Physics Working Group white
paper [323].

This section follows a similar structure to Sec. 4 where various modelling
approaches for compact binaries are discussed. However, the modelling of
gravitational waves emitted by cosmic strings is discussed in a separate
subsection 6.7, mainly because these sources require modelling approaches that
are specific to cosmic strings.

###  6.1 Numerical relativity

####  6.1.1 Beyond GR

NR simulations of BBH mergers in beyond-GR theories are only at their infancy
and have only been performed in a handful of theories, e.g. [1583, 638, 614,
1584, 1585, 630, 1586, 1587, 1588, 1589, 1590, 1591, 1592, 1593]. Even though
most studies have so far focused on proof-of-principle simulations, in recent
years there has been significant improvements in this topic.

The earliest NR simulations of BBH mergers in beyond GR theories focused on a
particular class of scalar-tensor theories that are known to possess a well-
posed initial value problem [1594, 1583, 638, 614]. In these scalar-tensor
theories, BHs satisfy no-hair theorems [392, 1595], and BBH mergers are
indistinguishable from GR. However, in more complex extensions of GR, or for
scalar matter with nontrivial boundary conditions (e.g., cosmological boundary
conditions or large scalar-field gradients), BHs can carry a scalar charge,
and these solutions need not be unique.  One of the main difficulties when
trying to do NR simulations in more complex theories is the fact that the
well-posedness of the initial value problem (see e.g. [509, 1596, 1597]) in
alternative theories is for most cases particularly challenging to establish
[1598, 1599, 1600, 1601, 1602, 1603, 1604, 1605, 1606, 1607, 1608, 1609, 1610,
1611, 1612, 1591, 1592, 1613, 1614].

To circumvent these problems, two treatments have been proposed in which
alternative theories are considered as effective field theories of gravity.
One approach uses an order-by-order expansion [1615, 1616, 1584, 630, 1586,
1587, 1617], in which the spacetime metric and extra fields are expanded
around GR, order-by-order in the coupling constant 444Sometimes this is also
called “order-reduction scheme.” This is not to be confused with order-
reduction schemes in which the field equations themselves are used to replace
higher-order curvature terms that are then kept only to a given (reduced)
order in the coupling constant., to guarantee well-posedness. The second
approach, the so-called equation-fixing method, proposes to use a scheme
inspired by the Israel-Stewart treatment of viscous relativistic
hydrodynamics, in which an effective damping controls higher frequency modes
while preserving the physics of the low frequency modes (nearly) untouched
[1601, 1618, 1619, 1620, 1593]. In this second approach the alternative theory
equations of motion are viewed as a low-energy effective deviation from the
Einstein equations, such that the theory can be captured by the low-energy
degrees of freedom. Although these methods are quite generic, they have so far
been mainly employed in theories with higher derivative terms in the metric
[1615, 1616, 1584, 630, 1586, 1587, 1621, 1619, 1612, 1617, 1622, 1593], such
as quadratic theories of gravity with an extra scalar field (for which the
order-by-order expansion was used), namely dynamical Chern-Simons gravity
(dCS) [404, 1623, 1624, 1625], Einstein dilaton Gauss-Bonnet gravity (EdGB)
[403, 1626, 1627, 334] or more generic Einstein scalar Gauss-Bonnet theories
(ESGB) [1628, 1629, 1630]. Unlike the simplest scalar-tensor theories, in dCS
and EdGB BHs naturally possess scalar “hair”, and therefore deviations from GR
naturally occur. Some ESGB theories are prone to BH scalarization, i.e., a
process in which GR BHs become unstable against the spontaneous development of
scalar hair for large enough spacetime curvatures [1628, 1629, 1630, 1631]
(see also [1632, 1633, 1634, 1635, 1636, 1637, 1622] for a similar process but
where scalarization is controlled by the magnitude of the BH spin).

In particular, using the order-by-order expansion mentioned above, BBH
simulations were done in dCS [1586, 1587] and EdGB [1621, 630]. Importantly,
this scheme also simplifies the problem of covering the parameter space needed
to build IMR waveforms for parameter estimation purposes. Since NR is scale-
invariant and the dCS and EdGB coupling parameters are scaled out of the
order-by-order expansion scheme [1586], a simulation for one set of BBH
parameters (masses and spins) can be applied to any total mass, and any valid
value of the coupling constants. With enough NR simulations, one can thus hope
to cover the the full BBH parameter space, build a surrogate model for
parameter estimation [590], and perform model-dependent tests of GR with
gravitational waves. Work towards model-dependent tests using full IMR signals
obtained through this method has been done in Ref. [1638], where parameter
estimation was performed on full BBH waveforms in dCS.

Despite their promises and advantages, the schemes mentioned above also suffer
from some important problems. The main problem with the order-by-order
expansion scheme is the fact that solutions may be plagued by secularly
growing solutions signaling the breakdown of perturbation theory on long
timescales. As proposed in Ref. [1639], such problem could potentially be
solved using renormalization-group techniques that can be used to build
solutions valid over the secular timescale. For the equation-fixing method on
the other hand, there are first studies in spherical symmetry [1619, 1640,
1620, 1641] and for BBH mergers in some beyond-GR theories [1593]. A more
detailed quantitative analysis is needed to identify how much a given
“equation-fixing” choice affects low-energy modes in scenarios of interest,
and to fully characterize the accuracy of these approximate schemes [1642].

These problems can be avoided for cases where a given theory in its full form
can be shown to have a well posed initial value problem. Important work in
this direction has recently been done. For example, nonlinear evolutions of
BBHs in quadratic gravity were presented in Ref. [1643]. Furthermore, BBHs
simulations have been performed in cubic Horndeski theories [1589] which are
known to be well-posed in the standard gauges used in numerical GR [1605]. For
more generic theories, it has been recently shown that, at sufficiently weak
coupling, the equations of motion for Horndeski gravity theories possess a
well posed initial value problem in a modified version of the generalized
harmonic gauge formulation [1606, 1607]. Quite remarkably, using this
formulation it could in principle be possible to perform simulations in all
Horndeski theories of gravity without having to resort to the two schemes
mentioned above [1592]. In particular, Ref. [1588] used this modified
generalized harmonic (MGH) formulation to perform the first numerical
simulations of BBH mergers in shift-symmetric ESGB gravity, without
approximations, which were also recently compared against PN results in Ref.
[1644]. Simulations of head-on BH collisions in ESGB theories that exhibit
spontaneous BH scalarization [1590] and simulations of binary neutron stars in
shift-symmetric Einstein-scalar-Gauss-Bonnet [1645] have also been done. Based
on the MGH formulation, a modified version of the CCZ4 formulation of
Einstein’s equations has also been proposed, which has been used to prove the
well-posedness of the most general scalar-tensor theory of gravity with up to
four derivatives, in singularity avoiding coordinates [1591]. This was then
used to perform (circular and equal-mass) BBH merger simulations in such
theories [1591, 1646].

Finally, NR simulations have also been done in theories in which BHs possess
an electric or magnetic charge [1647, 1648, 1585, 1649, 1650, 1651], possibly
coming from mini-charged dark matter [395]), primordial magnetic monopoles
[396]), or in specific classes of scalar-tensor-vector modified theories of
gravity [394]. In all these cases the field equations can be cast in a form
that is mathematically equivalent to Einstein-Maxwell theory (or to Einstein-
Maxwell with an extra scalar field for some theories [1585, 1652]), which is
known to possess a well-posed initial boundary value problem [1653].

####  6.1.2 Boson clouds

Simulations of BHs surrounded by massive bosonic fields have been widely
studied, focusing mostly on the evolution of bosonic fields around isolated
BHs (see e.g. [352, 1654, 1655, 366, 1656, 367, 1657, 1658, 356, 357, 364,
368, 373, 355, 1659, 1660]). This has led to the first successful nonlinear
evolution of the superradiant instability of minimally-coupled Proca fields
around a spinning BH in [364, 368]. These works have laid down the possibility
to study the impact of minimally-coupled bosonic fields in BBH systems within
NR [1661, 374, 378, 379, 1662, 375, 1663] but such studies are still in their
infancy.

####  6.1.3 Exotic compact objects

ECOs come in many different flavors, and the state of numerical simulations of
ECO mergers vary widely, depending on the source. For example, simulations of
string-theoretical models including fuzzballs, wormholes, firewalls or
gravastars have not yet been performed because the theoretical foundations —
the field content, equations of motion, equations of state (where applicable),
solution space — are under active investigation. Formulations suitable for
numerical relativity are not yet available.

In contrast, binary boson stars that are compact objects composed of real or
complex bosonic (scalar or vector) fields, with and without self interaction
potentials, have received a lot of attention [411, 1664, 1665, 1666, 1667,
1668, 438, 1669, 1670, 1671, 1672, 1673, 1674, 1675, 1676, 1677, 1678, 1679,
1680, 1681, 1682, 1683, 1684]. These boson stars may be treated as potential
dark matter components (see Sec. 2.7), or simply as toy model proxies for
unknown exotic objects. Several works have also considered mixed collisions of
BHs or neutron stars with these types of bosonic ECOs [1685, 1686, 1687].
Besides these simulations, mergers of extremely compact perfect fluid stars
have also been performed in [1688, 1689]. Most binary boson star simulations
have explored new phenomena that arise during their merger or the qualitative
structure of the gravitational radiation emitted. For model specific tests of
gravity more accurate waveforms will be needed which, in turn, requires
improved initial data and more accurate numerical evolutions.

In summary, despite the progress in the last few years, most of the studies
mentioned above have only provided a proof of principle for the stable
evolution of some families of ECOs and of fundamental fields around BHs in NR,
and have primarily aimed at identifying significant qualitative differences
between their signals and those of traditional binary candidates. Key areas
for improvements include the numerical accuracy of such simulations and the
quality of the initial data [1690, 1691, 375, 1683, 1663]. For example for
boson stars, a formalism is currently lacking to ensure that the superimposed
boosted objects are not in an excited state, the risk being that such
artifacts in the initial data are wrongly attributed as features of the ECO
signal. Less challenging, but equally important refinements are reducing the
eccentricity of the initial orbits (most simulations thus far have been head-
on mergers, or only achieve a few orbits before the plunge), and initial
errors in the Hamiltonian and momentum constraints (many simulations rely on
constraint damping which reduces control over the initial data parameters).
Efforts to solve some of these problems are currently underway, in particular
see [1691] where a novel method to solve the initial data constraints has been
formulated and that could be particularly useful when the sources involve
fundamental fields. See also Ref. [1684] where numerical simulations of binary
boson stars were performed using constraint satisfying initial data.

Whilst feasible given these improvements, a significant effort would be
required to create a template bank of waveforms of a similar quality to those
used for BH and neutron star mergers. This is true even if the class of
objects were restricted to a simple model, such as minimally coupled, non
self-interacting, non spinning complex scalar boson stars. Refining the NR
tools with which we study such objects is work that is merited in the run up
to the mission. In particular, key goals are (1) unambiguously identifying
characteristic deviations in the merger signals of ECOs and of BHs surrounded
by fundamental fields (e.g. to be used to motivate parameterised deviations in
IMR consistency tests) and (2) expanding the classes of objects for which
simulations of mergers have been performed to qualitative examples in all
sufficiently predictive ECO cases. A first NR waveform catalog of Proca star
mergers is available [1681].

###  6.2 Weak-field approximations

####  6.2.1 Beyond GR

Applying weak field approximation techniques to theories beyond GR raises
several difficulties which have prevented obtaining high-order PN results,
with the exception of a small set of alternative theories of gravity. For most
cases only the leading non-GR correction is known (see e.g. [400, 1692]). The
most advanced results concern scalar-tensor theories, for which the equations
of motion are known up to 3PN order [1693, 1694, 1695], including finite-size
effects [1696, 1697]. The gravitational flux and waveform are known up to 2PN
order [1698] while the scalar energy flux is known to 1.5PN order beyond the
quadrupolar formula [1699, 1700]. We expect the full scalar plus gravitational
waveform to be known up to 2PN order beyond the GR leading term in the next
few years, which should result in the first full waveform model beyond GR
[1701] shortly thereafter. The scalar-tensor equations of motion have been
used to obtain similar results in other theories such as Einstein-Maxwell-
scalar [1702, 1703] and ESGB theories [1704, 1705, 1706, 1697]. The main
difficulty encountered in more involved theories, such as the Horndeski family
[1707], is related to the need for a mechanism to screen the fifth force in
the solar system. As an example, the Vainshtein mechanism requires to keep
non-linearities in the description of the dynamics even in the weak field
limit, rendering the definition of a perturbative method challenging [1708,
1709, 1710, 1711, 1712, 1713, 1714, 1715].

Given the many different alternative theories of gravity proposed over the
years and the difficulty in constructing accurate waveforms in many of these
theories, the most common approach to build beyond GR waveforms for data
analysis purposes has been to use theory-agnostic phenomenological frameworks,
such as the parametrised post-Einsteinian (ppE) formalism [1716], the
generalized IMR formalism [1717] or the flexible theory-independent approach
[1718]. In these frameworks, deviations from GR are assumed to be small and
are treated as perturbative corrections to the signal predicted by GR. For
example, in the ppE formalism the inspiral phase of the frequency domain
waveform is given by:

$\tilde{h}(f)=\tilde{h}_{GR}(f)[1+\alpha_{i}(\pi\mathcal{M}f)^{a_{i}}]e^{i\beta_{j}(\pi\mathcal{M}f)^{b_{j}}}\,,$
(22)

where $\tilde{h}_{GR}(f)$ is the frequency domain waveform of GR,
$\mathcal{M}$ is the chirp mass of the binary, while the deviations from GR
are described by the dimensionless parameters $\alpha_{i}$,$\beta_{j}$,
$a_{i}$ and $b_{j}$. The index $i$ and $j$ indicate the PN order at which the
modification enters and are also introduced as a reminder that modifications
could enter at different PN orders, in which case we should sum over all of
them. Different values for the generic beyond GR parameters can be mapped to
distinct physical effects and gravitational theories (see e.g. [1719, 1720,
1721, 1722, 1723, 1724, 1725, 1726, 1575, 1727, 1728]), allowing to use this
formalism to place constraints on such theories. Extensions of the ppE
formalism to include higher-harmonic GW modes have also been recently
formulated [1729, 1718], as well as an extension for precessing binaries
[1730]. The latter was based on the recent computation of time- and frequency-
domain analytical waveforms emitted for quasicircular precessing BH binaries
in dCS gravity [1731] and in binaries composed of deformed compact objects
with generic mass quadrupole moments [469].

As a final remark, we should note that while the ppE approach is a very useful
and powerful method to perform null tests of GR, a qualitative and
quantitative interpretation of a constraint on the ppE parameters (or
measurements of a deviation from GR) only makes sense when accessed with
respect to specific theories [1732]; see e.g. Refs. [1733, 1734]. Hence
theory-agnostic approaches should only be seen as complementary, and not as
alternatives, to explicit computations of waveforms in specific alternative
theories.

####  6.2.2 Dark matter and boson clouds

Weak field approximations have been widely used to study the impact of dark
matter environments on binary systems [1735, 64, 1736]. Such an environment
can affect the binary dynamics through dynamical friction [1737], accretion
and the gravitational pull of the environment itself [1735, 64, 1736]. These
effects typically modify the waveform at negative PN orders and can be
captured within the ppE formalism presented above [64, 1736].

In this context, Newtonian approximations have also been used to study binary
systems moving inside bosonic structures, mainly motivated by the fact that
ultralight bosons are promising candidates to describe dark matter [337]. For
example, Refs. [1738, 1739] studied how compact binaries travelling through
these structures would be affected by dynamical friction and emission of
scalar radiation. In particular, it was found that, for sources in the LISA
band, scalar radiation affects the gravitational waveform at leading -6 PN
order with respect to the dominant quadrupolar term [1738, 1739]. In addition,
BHs surrounded by boson clouds have been studied within a weak-field
approximation [380, 381, 382, 1740, 372, 383, 374, 1741, 1742, 1743] showing
that a plethora of signatures, including tidally induced resonances [380,
382], floating and kicking orbits [1744, 1740, 382, 1745] and ionization of
the cloud [1742, 1746, 1747, 1748] can occur. Finally, it was recently shown
that sufficiently light bosonic fields can also be bound to and engulf a
binary BH system as a whole, showing a rich phenomenology [1749, 377, 378],
which could potentially lead to additional signatures besides the ones we
discussed above. A main problem with those studies is that they have mostly
been performed within a weak-field approximation, neglecting high PN
corrections and nonlinear effects close to merger. It would be especially
important to extend such studies beyond the weak-field regime and further
explore the physics of ultralight fields around binary BHs with full NR
simulations (see Sec. 6.1.2 for recent attempts).

####  6.2.3 Exotic compact objects

In the inspiral phase, the nature of the coalescing bodies can be studied
through (i) their multipolar structure; (ii) (the absence of) tidal heating;
and (iii) their tidal deformability. In particular, the multipole moments of
an ECO will be different from those of their Kerr counterpart [459, 460, 461,
462, 463, 464, 465, 466, 467, 468, 469, 470], which at leading order enters as
a 2PN correction to the GW phase due to the object’s quadrupole moment [1750,
469]. In addition, tidal interactions during the coalescence also provide
unique signatures able to disentangle BHs from ECOs in the form of tidal
heating. For BH binaries a small fraction of the emitted radiation is lost
through the horizon [1156, 830, 1751, 1752, 1753]. Absorption at the horizon
introduces a 2.5PN (4PN) $\times\log v$ correction to the GW phase of spinning
(nonspinning) binaries, relative to the leading term. On the other hand,
exotic matter is expected to weakly interact with gravitational waves, leading
to a smaller absorption during the inspiral, and therefore to a suppressed
contribution to the accumulated GW phase from tidal heating [449, 1754, 1755,
1756, 1757]. Finally, tidal deformations can be strong enough, especially
during the late stages of the inspiral, to modify the binary’s orbital
evolution leaving an imprint on the emitted waveform encoded in the ECO’s
tidal Love numbers [1758], which at leading order introduces a 5PN correction
to the GW phase [1036]. Overall more work is needed to construct fully
consistent inspiral waveforms for ECOs that incorporate all these ingredients.
For example, for the inspiral of boson stars with quartic interactions, such a
waveform was only recently constructed in Refs. [1759, 1760]. It would be
important to extend this to other types of ECOs.

###  6.3 Perturbation theory for post-merger waveforms

####  6.3.1 Beyond GR

The ringdown phase of a GW signal emitted by a BBH merger event can be
described using BH perturbation theory and is dominated by the QNMs of the
remnant [1327, 1329]. Owing to the very large SNR with which LISA will be able
to detect the ringdown phase of supermassive BHs [612, 1358], we expect to be
able to measure several QNMs for a single event up to very large distances
[1335, 612, 1338, 1761]. This will allow to perform precise consistency tests
of the QNMs and test the hypothesis that the remnant of a BBH merger is well
described by the Kerr metric in GR [1335, 612, 1358, 1762, 1763, 1764, 1765,
1766, 1369, 1336, 1767, 1374, 1768, 1769, 1367, 1371, 606, 1770, 1364, 1368].
Due to these prospects, significant work has been done in recent years in
order to accurately model the ringdown phase of BBH mergers within GR [1438,
1361, 1363, 1771, 1369, 1370, 1772, 1773, 1354, 1373, 1368, 1376, 1393].

For beyond GR theories, however, there has been significantly less progress,
even though QNM measurements can, in principle, be powerful probes of beyond
GR physics. In particular, for beyond GR theories, the ringdown of a BBH
merger is typically expected to differ from the ringdown in GR (although not
necessarily for all non-GR theories). Computations of QNMs in beyond GR
theories are very challenging, especially for rotating BHs for which the
separability of the perturbative equations present in GR [1324] is in general
absent. Then, one may employ numerical methods to compute QNMs for non-
separable equations, see e.g. Refs. [1774, 1775, 1776, 1777, 1778]). In
addition, most rotating BH solutions in alternative theories, for the cases
where they differ from Kerr, are only known either perturbatively, through a
small-spin expansion around non-spinning backgrounds [1779, 1780, 1781, 1782,
1783], or they are given in the form of numerical solutions of the field
equations (e.g. [1784, 359, 1785, 360, 1786, 1787, 1788]), complicating the
problem even further. Furthermore, BH spacetimes may change their character as
compared to GR. For example, it was shown that BHs in ESGB and dCS gravity are
of Petrov-type I [1789] (instead of the more symmetric type-D spacetimes like
the Kerr metric). Therefore, even though QNMs have been computed for a handful
of theories and BH solutions (see e.g. [1790, 1791, 64, 1792, 459, 1793, 1794,
1782, 1795, 1796, 1797, 1798, 1799, 1800]), the vast majority of these results
either assumed non-rotating or slowly-spinning BH backgrounds or relied on
approximations such as the eikonal/geometric optics approximation. Methods to
circumvent some of these problems were formulated in Refs. [1801, 1777, 1802,
1803, 1804, 1805] where modified Teukolsky equations governing the
perturbations of non-Kerr spinning BHs arising in theories beyond GR, were
derived in the case where the deviations from the Kerr geometry are small.
Another approach would be to rely on fits of time-domain waveforms computed
from numerical simulations in selected classes of theories beyond GR. However,
this is quite a challenging task, even in GR [1361, 1362, 1363, 1771, 1369,
1370, 1772, 1773, 1375, 1354, 1373, 1376, 1393], besides the fact that such
simulations are still only possible for a handful of cases, has emphasized in
Sec. 6.1.1.

####  6.3.2 Exotic compact objects

The vibration spectra of ECOs have been computed for a wide class of models,
although mostly for spherically symmetric configurations, including: boson
stars [1806, 1807], gravastars [1808, 1809, 1810, 1811], wormholes [437, 1812,
1813], quantum structures [1814, 1815, 438, 1816, 1817, 1818, 1819], and in a
model-indepedent fashion using the membrane paradigm [436]. Typically, the
QNMs of ECOs differ from the BH QNMs due to the presence of a surface instead
of an event horizon [102] and the excitation of the internal oscillation modes
of the objects [1820, 1821, 459]. In addition, the isospectrality of axial and
polar modes of spherically symmetric BHs in GR [1326] is broken in ECOs, which
are instead expected to emit a characteristic “mode doublet”. The detection of
such doublets would be a irrevocable signature of new physics [1329, 1817,
1797, 436]. The formation of an exotic ultracompact object may also lead to
the emission of GW echoes in the post-merger signal [437, 434] (see also Sec.
2.7). Several phenomenological GW templates for echoes have been developed
based on standard GR ringdown templates with additional parameters [1822,
1823, 433], the superposition of sine-Gaussians [1824], hybrid methods that
put together information from perturbation theory and the pre-merger orbital
dynamics [1825, 1826, 1827], using the close-limit approximation [1828], and
analytical models that explicitly depend on the physical parameters of the ECO
[1829, 1830, 456, 1831].

###  6.4 Small mass-ratio approximation

####  6.4.1 Beyond GR

Although there have been substantial advances in GSF models in GR (see Sec.
4.3), self-force calculations for theories beyond GR are in their infancy. The
most detailed studies in this vein have so far investigated changes to the
emitted GW flux for scalar-tensor and higher-order curvature theories of
gravity [134, 1832, 1744, 1833, 1834, 1792, 1835], and for theories of massive
gravity when assuming a Schwarzschild BH background [1836]; there are only a
few examples of formulations of a full self-force theory beyond GR [1837,
1838]. But recently it has been argued that for a vast class of theories, no-
hair theorems or separation of scales lead metric and scalar perturbations to
decouple [1839]. In particular, in the large class of higher-derivative
gravity models, the modification to GR scales with the curvature, i.e., with
the inverse of the mass. Consider, for example, the Kretschmann curvature
scalar $\mathcal{K}=R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$. On the horizon
of a Schwarzschild BH of mass $M$ it scales as $\mathcal{K}\sim M^{-4}$.
Consequently, higher-curvature modifications to a supermassive BH are
negligible. This allows the background spacetime to be treated as the Kerr
solution, and changes in the binary’s evolution to be fully controlled by the
scalar field charge of the small BH [1839, 1840, 1841, 1842, 1843]. In this
regard, the study of small-mass-ratio binaries in generalized scalar-tensor
theories of gravity can benefit from efforts already devoted to investigating
self-forces on scalar charges, and their orbital evolution around BHs [1844,
1845, 1846, 1847, 1301, 1848, 1849, 1850, 1851, 1852, 1240, 1069, 1098, 1241].

####  6.4.2 Dark matter and boson clouds

For EMRI and IMRI systems, the largest effect of dark matter is to produce a
significant drag force through dynamical friction [1737], which increases the
rate of inspiral of the small compact object [376, 64, 1853, 1736, 1854]
(accretion and the usual gravitational pull of the dark matter can be
significant, but are in general subleading effects [1735, 64, 1855, 1736]).
This typically requires very high densities of dark matter surrounding the
intermediate or supermassive BH [64, 1853]. The required densities can be
achieved when the IMBH or SMBH adiabatically grows in a dark-matter halo, and
the distribution of dark matter near the BH gets compressed into a “spike”
with significantly higher densities than the surrounding halo [1856, 1857,
1858, 1859]. Sufficiently large densities may also be achieved for clouds of
ultralight bosons formed around MBHs, as discussed in Sec. 2.7. If they exist,
these clouds would introduce additional metric perturbations that perturb the
companion’s orbit, and they would respond to the gravity of the companion. If
the cloud contributes a sizable fraction of the BH’s mass, its gravity would
necessitate a change in the background geometry (see Refs. [1860, 1861] for
such an example), dramatically complicating the GSF model. Building accurate
IMRI/EMRI waveforms evolving in a dark matter spike or boson cloud is also
crucial to distinguish these environments from other astrophysical
environments, such as accretion disks. In fact, recent studies strongly
suggest that these different types of environments leave distinguishable
signatures in IMRI or EMRI waveforms [384, 385, 1862, 1289, 1863, 1864, 1859].
Therefore the detection of such systems could also be used to identify and
learn about the environment surrounding supermassive BHs.

To date, studies of the gravitational effects of dark matter spikes and boson
clouds on binary systems, such as EMRIs and IMRIs, have mostly relied on
Newtonian or post-Newtonian approximations for the orbit and the matter
distribution [1735, 1853, 1855, 1865, 384, 385, 1745, 383, 380, 1745, 1740,
372, 1866, 1746, 1747, 1867, 1859, 1748]. Generalizing some of these
calculations to the fully relativistic case is in principle feasible. Indeed,
first steps in this direction were recently done in Refs. [1868, 1291, 1290]
where a generic and fully-relativistic formalism to study EMRI systems in
spherically-symmetric, non-vacuum BH spacetimes was developed. For boson
clouds, a fully relativistic framework that assumes that the impact of the
cloud on the BH geometry can be treated perturbatively, was proposed in Ref.
[1869].

However, there are important open problems to be solved. First, generalizing
the formulation of Refs. [1868, 1291, 1290, 1869] to non-spherically-symmetric
backgrounds is a non-trivial task, and secondly, for dynamic matter
distributions [1866], evolving the matter environment coupled to the binary
system at first post-adiabatic order will likely be a complex problem.

Finally, it is worth noting that boson clouds themselves can also be strong
sources of nearly-monochromatic gravitational waves potentially detectable by
LISA [335, 332, 365, 342, 387, 343, 389, 369, 371, 370, 1870]. The
gravitational waves emitted by these sources can be computed using the same
techniques typically used in the small mass-ratio approximation: one can
consider the cloud as being a small perturber of a Kerr background geometry
and use the Teukolsky formalism to compute GW emission from a boson cloud
[365, 343, 369, 371, 1870].

####  6.4.3 Beyond GR BHs and exotic compact objects

A key goal of the LISA mission, is to determine whether the dark central
objects in galactic cores are genuine BHs, and if so, whether they are
accurately described by GR BHs. With EMRIs we expect to be able to constrain
fractional deviations from the quadrupole moment of the Kerr solution at a
level smaller than $10^{-4}$ [105, 99], which would allow to impose stringent
constraints on ECOs and non-GR BHs. Work to include these modifications in GSF
models has been done by considering the large central object as a “bumpy” BH.
Bumpy BHs are spacetimes that include the Kerr limit, but that differ in a
parameterized way, typically by modifying the spacetime’s multipole moment
structure [1871, 104, 1872, 1873, 1874, 468]. The “bumps” can be introduced
directly into the background metric, or (if they are sufficiently small) they
can be treated as additional metric perturbations. A benefit of the bumpy BH
framework is that it is agnostic about the mechanism which produces the BH’s
bumps. They could describe non-GR physics (to exactly recover BH solutions in
specific alternative theories [1875, 1876], for example), or if the central
object is an ECO, they could be tidally induced by the companion, encoding the
ECO’s tidal Love numbers [1877, 1282, 1878, 1879, 1880, 1881, 1758, 1882, 471,
1883, 472]. Because bumpy BHs break the symmetries of Kerr that yield
integrable motion [1884], orbits in these spacetimes are generically chaotic
[1885, 1886, 1887, 1888] and subject to prolonged resonances [1889, 1887,
1318, 1888]. Most of the work considering EMRIs around spacetimes with a
multipolar structure different from Kerr, constructed approximate “kludge”
waveforms [1890, 488, 1874, 1891, 468] generated considering geodesics in a
perturbed Kerr spacetime with parameters evolved using post-Newtonian
equations [107, 105, 1885, 1233, 1873]. While these approaches are believed to
reproduce the main features of the orbit, they will not be enough to get to
the precision required to determine that the inspiral is indeed an inspiral
into a Kerr BH or not [1081, 1885], and much more work is needed to build such
waveforms. The additional parameters in such models will also typically
increase the possibility of degeneracies, making it important to fully
understand possible discernible features, such as characteristic variations of
the amplitude and the energy emission rate [1264] or the appearance of
prolonged resonances [1233, 1318, 1888]. If the central object is an ECO, in
addition to nonintegrable motion and geodesic resonances [1892], it will also
lack a true event horizon. This causes a change in boundary conditions: in a
BH spacetime, fields are purely ingoing on the horizon, while in an ECO
spacetime, fields obey a boundary condition of (partial) reflection at the
effective surface. The effect of this reflection has been taken into account
for the dissipative first-order self-force in [1893, 1894, 1756, 1895, 457,
1896]. In addition, the modes of the ECO could also be excited during the
inspiral of the small compact object [1897, 1807, 1894, 1898, 1899, 1895, 457,
1900] (see Ref. [1900] where the conditions for such modes to be effectively
excited during the inspiral were studied in more detail). Depending on the
object’s properties, this could cause the rate of the inspiral to change
significantly in the vicinity of each resonance. For the dissipative first-
order self-force, the effect on the emitted GW flux has been modelled for
gravastars [1897], boson stars [1807] and an exotic ultracompact object with
an exterior Schwarzschild or Kerr geometry but with a hard surface close to
the would-be horizon [1894, 1895, 457]. An object plunging in the interior of
an ECO would also be subject to additional forces that can dominate the
dynamics, such as dynamical friction [1737] and accretion. The impact of these
effects on the orbit dynamics and GW waveforms has been modeled for small
compact objects plunging onto boson stars using a Newtonian approximation
[376] (see also [1901] for earlier work on the same subject were dynamical
friction and accretion had been neglected).

###  6.5 Effective-one-body waveform models

EOB waveforms with generic deviations from GR have been developed based on a
parameterized form for performing theory-agnostic tests. These include
deformations to the GW phase (see Eq. (22)) which are added as corrections to
frequency-domain EOB waveforms in GR [1902, 1903, 1718, 604, 605], or
parameterized deformations of the QNMs in time-domain EOB waveforms, which
have been used to perform tests of the no-hair theorem in GR and constrain
higher-order curvature theories of gravity [1764, 606, 1904, 604, 605, 1905].

Recent work has also computed the EOB dynamics for non-spinning binaries in
various theories beyond GR as a foundational step towards enabling model-
dependent tests. Results are available for scalar-tensor theories of gravity
[1514, 1515, 1516, 1906], Einstein-Maxwell-dilaton theories [1907] and
Einstein-scalar-Gauss-Bonnet gravity [1704, 1906]. The EOB approach has also
been used to compute the two-body potential energy of point-like particles in
theories exhibiting a Vainshtein screening [1908]. Developing these results
into full EOB waveform models and including more realistic physics remains an
open task.

For non-vacuum binaries involving e.g. ECOs or boson clouds around BHs, a
number of generic matter effects during the inspiral change the GW signals
away from those of a BH binary (see discussion in previous subsections).
Several of the dominant matter effects for neutron star binaries are included
in EOB models that could also have broader applicability to ECOs or non-vacuum
BHs. These effects are described in a parameterized form, where the nature of
the object is encoded in characteristic coefficients such as the bodies’
multipole moments, tidal and rotational Love numbers, and QNM frequencies.
Specifically, the effects of rotational deformation and adiabatic tides are
currently included in the TEOBResum family of models discussed in Sec. 4.5 in
the way described in Refs. [1909, 1038, 1910, 1911, 1446, 1442, 1452, 1481,
1475, 1484], with the state-of-the-art models incorporating the known
adiabatic gravitoelectric and -magnetic tidal terms up to $\ell=8$ in
different resummation schemes [1910, 1911, 1912, 1481], dynamical tides from
the fundamental modes [1913, 1480, 1479] as well as nonlinear-in-spin effects
dependent on the nature and interior structure of the object [1442, 1475]. The
SEOBNR models include tidal effects using the dynamical tidal models for the
fundamental modes from [1480, 1479], which have recently been further
developed [1482, 1914, 1915]. Matter effects have also been incorporated in
the ROMs of SEOBNR waveforms in the frequency domain by augmenting the ROMs
developed for BBHs with analytical closed-form expressions correcting the GW
phase to include tidal and spin-induced multipole effects based on PN
calculations for the rotational quadrupole effects [1916] and tidal models
that are specific to calibrations to NR simulations of neutron star binaries
[1917, 1533], though in principle could be replaced by more general frequency-
domain tidal models.

However, significant work remains. Besides the improvement of current EOB
models already discussed in Sec. 4.5, other improvements that are specific to
ECOs or non-vacuum BHs include, for example, incorporating more effects from
spins and relativistic phenomena in the matter sector (which are likely to be
more important for compact ECOs than for neutron stars); incorporating
parameterized absorption coefficients; and the inclusion of effects specific
to dark matter spikes and boson clouds already discussed in Secs. 6.2.2 and
6.4.2. In addition, a remaining challenge is to develop full IMR waveforms.
Possible avenues for constructing full waveforms for binaries with matter have
been explored, e.g., in the context of binary neutron stars [1918], by using
NR-informed analytical models that build on quasiuniversal features found in
these systems [1919, 1918, 1920, 1920, 1920], and for the tidal disruption in
BH-neutron-star binaries [1921, 1483, 1484]. To achieve similar complete
waveforms for ECOs, it is worth mentioning that a phenomenological model was
proposed in Ref. [1922], which used EOB waveforms with tidal effects to
describe the inspiral phase and a toy model inspired by Ref. [1923] to model
the post-merger dynamics of the system. Whether this model is an accurate
enough description for known ECOs, e.g. boson stars, remains to be fully
studied.

###  6.6 Phenomenological waveform models

The phenomenological waveform models of Sec. 4.6 already provide standard
tools for tests of GR with the current generation of ground-based detectors
[1924, 1902, 1925, 604, 605]. In particular, parameterized tests (in the
spirit of ppE formalism, see Sec. 6.2.1) have been employed to probe theory-
agnostic deviations from GR in the Phenom ansatz describing the phase
evolution of an observed signal [1717, 1718] and to test the BH nature of
compact binaries from waveform signatures of spin-induced quadrupole moments
[1926]. At a more fundamental level, such tests depend on the parameterization
of the waveform model, and face a number of practical and conceptual problems,
e.g., varying the phenomenological coefficients in such models by amounts that
are too large, pathologies may arise in the waveform. Some improvements can be
expected from adding time domain models [1502] to the existing frequency
domain models, thus providing very different waveform parameterizations to be
varied. Going beyond theory agnostic tests, the flexibility of
phenomenological models also provides a natural ground to develop models for
beyond GR theories or “exotic” physics, such as boson stars and boson clouds
around BHs. A key problem is the very large parameter space for such theories,
and significant technical and conceptual progress will be required to develop
mature approaches. Natural first steps for the near future would be, for
example, to calibrate accurate models to small “toy model” regions of some
beyond GR theories, and to develop simple toy models that do not require any
calibration to numerical waveforms.

###  6.7 Modeling cosmic strings

Waveform models for GW emission are needed for LISA for three types of
anticipated searches: for the stochastic background, for individual bursts and
for individual, coherent, nonlinear oscillatory waveforms. There are
significant uncertainties in the properties of the underlying string elements
(the number of objects, the presence of cusps, kinks and small scale
structure, the size distribution of the objects, the spatial and velocity
distribution of the objects). Here, however, we concentrate on the modeling of
the waveforms generated by one or more loops while also mentioning some of the
incompletely known loop properties that will influence the waveforms.

####  6.7.1 Nambu–Goto, smooth loops with cusp/kink features

We begin with minimally coupled string networks with loops that are smooth on
large scales and include cusps or kinks but no other small scale structure.

The underlying formalism to calculate the gravitational signal to lowest order
in the string tension $\mu_{S}$ is well known [1927]. Early research focused
on calculating the individual loop waveforms for relatively simple loops
[1928, 1929, 282, 283, 1930, 1931, 1932, 1933, 1934]. Both the burst waveforms
and the relevant statistical averages for the background implicitly follow.
FFT techniques work well to estimate the power at low to moderate harmonics
[1930, 263] and additional numerical techniques have been devised specifically
to handle higher harmonics [1935].

In the case of LVK, small scale features dominate emission at high
frequencies. A similar situation pertains to LISA, which will be sensitive to
the high frequencies of the astrophysically relevant low tension strings. The
development of an asymptotic approach Refs. [269, 270, 298] was an important
quantitative and qualitative advance. High frequency emission sourced by cusps
and kinks is approximately independent of the large scale loop structure. The
important features of cusp waveforms in this limit are that (1) the cusp power
spectrum falls off at high frequencies as $f^{-4/3}$; and (2) the GWs from the
cusp are strongly beamed in the direction of cusp motion and exponentially
suppressed at frequencies $f\gtrsim 1/(\theta^{3}T)$ for angle $\theta$ with
respect to that direction (with $T$ the periodicity of the string). The
waveform’s time dependence is $\propto|t-t_{c}|^{1/3}$, where $t_{c}$ is the
time at which the cusp formation event is noted by an observer situated in the
beam of emission ($\theta=0$), and so the cusp waveform is called _sharp_ or
_spiky_. This sharpness is softened for $0<\theta\ll 1$. Similar analyses are
applicable to kinks.

The genericity of the high frequency predictions for cusps and kinks, the
independence with respect to the large scale loop structure and the simplicity
of the templates are all important, and the results have been exploited in
astrophysically-motivated searches [299, 1936, 302, 305].

The outlook and need for waveforms is as follows:

(1) The existing asymptotic treatments of kink and cusp emission are likely to
be sufficient for modeling an unresolved stochastic wave background. That is
because the energy flux derived from averaging over time, space and angle,
which is implicit in such a calculation, depends upon low-order moments of the
loop’s full beam pattern that can be derived from the asymptotic results.

(2) More refined models of the burst waveform that do not, strictly speaking,
lie in the asymptotic regime are needed. Such waveforms can account for
geometric effects when the observer does not lie directly in the cusp
direction and for anisotropies in the beam shape. These can be found by
extending the asymptotic treatment.

(3) Existing calculational methodologies for full loops are accurate and
efficient at low to moderate modes. LISA and pulsar timing arrays will probe
both low and high order modes. Hybridization methods giving models of the beam
that span the full range of modes will be required for any loop sources that
stand out above the stochastic background [1935].

####  6.7.2 Nambu–Goto, small-scale structure

Each time a long horizon-crossing string intercommutes it gains two kinks.
Small scale structure builds up on the string [1937, 1938, 1939, 1940, 1941].
Gravitational backreaction probably smoothens the string at small scales
[1942, 1943, 1944, 1945]. The wigglyness that remains has been evaluated
[1937, 1938] but never determined directly by simulation. Ref. [1946] found
that small-scale structure near the cusp rounded off the sharpness of the cusp
waveform but affected only observers very near to the direction of cusp
motion. Alternatively, Ref. [1938] found that the small scale structure
enhanced intercommutation on a newly formed loop and excised the cusp when it
first begins to form. It will be important to investigate the effect of small
structure on loop dynamics and the resultant waveforms. This can be done
parametrically as the amount and character of the wigglyness increases on
otherwise simple loops.

####  6.7.3 Nambu–Goto for Pseudocusps

As a last note, in addition to cusp waveforms, some authors [1947] have
studied the bursts due to “pseudocusps”, i.e., string trajectories which very