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

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

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

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

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

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

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

Coordinators: Marta Colleoni and Adam Pound

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

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

The perturbative field equations and trajectory of the secondary are determined using GSF theory, a collection of techniques for incorporating a small gravitating object into an external field [98, 1043–1046]. Using matched asymptotic expansions (or EFT [1047]), the small object is reduced to a point particle, or more generally to a puncture in the spacetime geometry, equipped with the object’s multipole moments. For a compact object, higher moments scale with higher powers of the small mass $m_2$, such that one additional multipole moment appears at each order in $\epsilon$. This multipolar particle (or puncture) is found to obey a generalized equivalence principle, behaving as a test body in a certain effective, smooth, vacuum metric $g_{\alpha\beta}^{(0)} + \epsilon h_{\alpha\beta}^{(1)} + \epsilon^2 h_{\alpha\beta}^{(2)} + \dots$ [1044, 1048–1060]. The regular fields $h_{\alpha\beta}^{(n)}$ that influence the secondary’s motion can be calculated directly, using a puncture scheme in which they are the numerical variables [1061, 1062]. Alternatively, they can be extracted from the full fields $h_{\alpha\beta}^{(n)}$ using a mode-by-mode subtraction [1063–1069] or other methods [1070, 1071], as reviewed in [1046, 1072, 1073]. In either approach, the perturbative field equations can be solved using the methods of BHPT [740, 1046, 1074–1083].

A key feature of this model is its clean separation of time scales: because the self-force is small, inspirals occur slowly, over $\sim 1/\epsilon$ orbits. On the time scale of a few orbits, the secondary’s trajectory is approximately a bound geodesic of $g_{\alpha\beta}^{(0)}$ (illustrated in Fig. 5), which in Kerr spacetime is generically triperiodic, undergoing radial, polar, and azimuthal motion with frequencies $\Omega_r$, $\Omega_\theta$, and $\Omega_\phi$ [1046, 1084–1087]. Over the long inspiral, the frequencies slowly evolve due to dissipation. The
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Geodesics of Kerr spacetime. Left: a non-resonant geodesic, which fills a toroidal shape. Right: a resonant geodesic. Images taken from [98].</figcaption></figure>

field equations can therefore be solved using a two-timescale expansion [1046, 1088–1093], $h_{\alpha\beta}^{(n)} = \sum_{k^A \in \mathbb{Z}} h_{\alpha\beta}^{(n,k^A)}(J^B, x^a)e^{ik^A\psi_A}$. Here $x^a$ are spatial coordinates, and all time dependence is encoded in the set of system parameters $J^A$ (the secondary’s orbital energy, angular momentum, and Carter constant, the primary’s mass and spin, etc.), and the set of orbital phases $\psi_A$. The parameters $J^A$ evolve slowly, on the inspiral time scale $t \sim 2\pi/(\epsilon\Omega)$, while the phases $\psi_A$ evolve on the orbital time scale $t \sim 2\pi/\Omega$, with evolution equations of the form [747, 1046, 1089, 1092, 1093]

$$ \frac{dJ^A}{dt} = \epsilon G_{(1)}^A(J^B) + \epsilon^2 G_{(2)}^A(J^B) + O(\epsilon^3), \quad \frac{d\psi_A}{dt} = \Omega_A^{(0)}(J^B) + \epsilon\Omega_A^{(1)}(J^B) + O(\epsilon^2). \quad (13) $$

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

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

Many of the core tools in GSF modelling, as well as advanced codes described in the sections below, have been consolidated in the open-source Black Hole Perturbation Toolkit [740], which provides a hub for GSF code development, as well as in the Black Hole Perturbation Club [1109]. In particular, the Fast EMRI Waveforms (FEW) package [1110], which exploits GSF models’ multiscale structure, has provided a flexible framework for rapid waveform generation.
**4.3.2. Suitable for what sources?** Historically, the GSF approximation has been motivated by EMRIs, with mass ratios in the interval $10^{-7} \lesssim \epsilon \lesssim 10^{-3}$ as far as LISA sources are concerned (see Table 4). The GSF model's accuracy for a typical EMRI can be estimated from the evolution equations (13). On the inspiral time scale, the solutions for the orbital phases, and therefore the GW phase, take the form

$$ \psi_A = \frac{1}{\epsilon} \left[ \psi_A^{(0)} + \epsilon \psi_A^{(1)} + O(\epsilon^2) \right]. \qquad (14) $$

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

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

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

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

**4.3.3.1. Adiabatic approximation** At leading order in the two-timescale expansion, only the geodesic frequencies $\Omega_A^{(0)}$ and time-averaged dissipative piece $G_A^{(1)}$ of the first-order self-force are required to drive the evolution [1085, 1123]. The secondary's motion in this approximation can be regarded as an adiabatic inspiral through a sequence of geodesic orbits [1089, 1124]. The waveform can then be built from a corresponding sequence of “snapshot” Fourier mode amplitudes together with the leading phases in
<table><thead><tr><th rowspan="2">Background Spacetime</th><th rowspan="2">Orbital Configuration</th><th colspan="2">Adiabatic</th><th colspan="3">Post-1-adiabatic</th></tr><tr><th>1SF (Dissipative)</th><th>1SF (Conservative)</th><th>2SF (Dissipative)</th><th>Spin Effects (Conservative)</th><th>Spin Effects (Dissipative)</th></tr></thead><tbody><tr><td rowspan="2">Schwarzschild</td><td>Circular</td><td>✓✓✓</td><td>✓✓✓</td><td>✓✓✓</td><td>✓✓✓</td><td>✓✓✓</td></tr><tr><td>Eccentric</td><td>✓✓✓</td><td>✓✓✓</td><td>X</td><td>✓✓, ✓✓✓*</td><td>✓, ✓✓*</td></tr><tr><td rowspan="4">Kerr</td><td>Circular</td><td>✓✓✓</td><td>✓✓</td><td>X</td><td>✓,✓✓*</td><td>✓✓✓*</td></tr><tr><td>Eccentric Equatorial</td><td>✓✓✓</td><td>✓✓</td><td>X</td><td>✓,✓✓*</td><td>✓✓*</td></tr><tr><td>Generic</td><td>✓✓✓</td><td>✓</td><td>X</td><td>✓</td><td>✓*</td></tr><tr><td>Resonances</td><td>✓✓✓</td><td>✓</td><td>X</td><td>X</td><td>X</td></tr></tbody><tfoot><tr><td colspan="2">✓✓✓ Evolving Waveform</td><td>✓✓ Driven Inspiral</td><td>✓ Snapshot Calculation</td><td colspan="3">*(Anti-)Aligned Spin Only</td></tr></tfoot></table>

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

Eq. (14). Alternatively, outside the two-timescale scheme, a time-domain field equation can be solved with an adiabatically inspiraling source particle [1117, 1125–1128].

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

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

Once fluxes and waveform mode amplitudes have been computed across the parameter space, adiabatic waveforms can be efficiently generated using a combination of neural network and reduced-order techniques [1110, 1160], a method that lends itself to graphical processing unit (GPU) acceleration techniques; see Sec. 5.2. This approach should ultimately meet the efficiency requirements described in Sec. 3.3.1. However, due to the high-dimensional parameter space, significant work remains to populate the parameter space. As a consequence, data analysis development has so far relied on semi-relativistic “kludge” models that approximate the adiabatic evolution and waveform generation [488, 1164–1167].
4.3.3.2. *Post-1-adiabatic self-force effects: 1SF* The 1PA terms in the two-timescale expansion ($\Omega_A^{(1)}$ and $G_{(2)}^A$) take as input the conservative piece of the first-order self-force (1SF, $\propto h_{\alpha\beta}^{(1)}$) and the dissipative piece of the second-order self-force (2SF, $\propto h_{\alpha\beta}^{(2)}$). Once these ingredients have been computed, a 1PA waveform-generation scheme takes the same form as a 0PA one [1046, 1092].

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

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

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

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

4.3.3.3. *Post-1-adiabatic self-force effects: 2SF* The dissipative piece of the second-order self-force ($G_{(2)}^A$) contributes to the GW phasing at the same 1PA order as the first-order conservative self-force, but calculations of it are less mature. After years of development of the governing formalism [1056–1060, 1220–1222] and of practical implementation methods [1046, 1073, 1091, 1223–1228], Refs. [1018, 1118, 1119] recently carried out the first concrete calculations of physical second-order quantities in the restricted case of quasicircular orbits around a Schwarzschild BH. These calculations culminated in the first complete 1PA waveforms in Ref. [1119]. Figure 7 shows a comparison between one of these 1PA waveforms and an NR waveform for $q = 10$. As alluded to in Sec. 4.3.2, the 1PA approximation agrees well with NR even at this moderate mass ratio.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 1PA GSF waveform for a quasicircular, nonspinning binary with mass ratio $q = 10$ (orange). The inset shows a zoomed portion of the waveform near the merger. Also included for comparison are the 0PA GSF waveform (blue, inset only) and the waveform for the same binary produced using an NR simulation in the SXS catalog (SXS:BBH:1107, in black). The three waveforms are aligned in time and phase at $t = 320M$, when the orbital separation is ≈ 13.83$M$. Image reproduced from Ref. [1119].</figcaption></figure>

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

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

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

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

### 4.3.3.6. Spin and finite-size effects

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

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

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

The secondary's higher moments are unlikely to contribute at 1PA order. The tidal quadrupole's contribution to the acceleration scales with $\epsilon^4$ [1282] and can only become relevant for non-compact objects [1267]. The spin also induces a quadrupole deformation of the secondary, which creates an acceleration $\sim \epsilon^2$ [87, 1122] that is unlikely important for EMRs (up to resonances) but might be relevant for IMRs [87].
4.3.4. *Environmental effects* In the preceding sections, we have focused on the simplest GSF model: an isolated vacuum BH orbited by a single body. This will not fully describe astrophysical small-mass-ratio binaries. Accreting matter will be present at some level, and there are likely to be additional bodies nearby. Accounting for these effects will be necessary to achieve the stringent accuracy requirements that EMRIs, in particular, place on GSF models.

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

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

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

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

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

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

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

Eccentricity and inclination bring multiple challenges. Most recent self-force calculations have utilized decompositions into Fourier modes (the natural setting of the two-timescale expansion) and angular harmonics. Eccentric, inclined orbits can require ~10⁵ modes, all of which will couple to one another in the second-order source. The sum of Fourier modes can also suffer from poor convergence: in the existing 2SF puncture scheme, the source for $h_{\alpha\beta}^{R(2)}$ has finite differentiability on the puncture's worldline, leading to slow power-law convergence (the Gibbs phenomenon). At first order, similar problems were overcome using methods of extended solutions [1171, 1172, 1177, 1181, 1300–1302] that restore exponential convergence. These have inspired a new scheme, applicable at second order, known as the method of extended effective sources [1303], which has been demonstrated in the case of a scalar-field toy model for eccentric orbits about a Schwarzschild BH. Work now remains to apply it to gravitational perturbations, both at first and second order, and to orbits in Kerr spacetime.
The extension to Kerr spacetime brings more challenges. Second-order calculations in Schwarzschild have so far relied on directly solving the perturbative Einstein equations, which are not separable in Kerr spacetime. At first order, this problem was overcome using radiation-gauge methods, in which the metric perturbation is reconstructed from a solution to the (fully separable) Teukolsky equation [1076, 1078, 1080]. One path to 2SF calculations is to extend that method to second order [1182, 1304–1306]. The standard method of metric reconstruction [1078, 1080] fails beyond linear order, but recent work found an extension to all orders [1182, 1305]. Recent work on second-order flux-balance laws [1307] also suggests that 1PA rates of change of energy and angular momentum (though perhaps not the Carter constant) can be computed directly from a solution to the second-order Teukolsky equation, without the need for metric reconstruction.

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

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

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

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

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

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

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

In addition to including more of the two-body GR parameter space, considerable work
must be done to include possible beyond-GR and environmental effects in GSF models.
Fortunately, the GSF model is relatively modular. This means most additional effects
can be added separately, as described in Sec. 6.4, and are readily incorporated into
frameworks such as FEW [1110]. However, modelling these additional effects will be
particularly challenging if they are not amenable to a two-timescale treatment, or if
they are too large to be treated perturbatively.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Real and imaginary parts of the Kerr quasinormal mode frequencies, coloured by the BH spin (from $a = 0$ in dark blue to $a = 0.9$ in yellow). Left: $\ell = m = 2$ fundamental and first overtone modes. Right: modes with different angular numbers $(\ell, m)$.</figcaption></figure>

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

Coordinators: Stephen R. Green and Laura Sberna

Contributors: E. Berti, R. P. Macedo, P. Mourier, N. Oshita, M. van de Meent

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

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

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

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

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

One can easily estimate the mass range of remnants with ringdowns within the LISA band. Indeed, the longest lived (fundamental, $\ell = m = 2, n = 0$) QNM has frequency $f_0 \approx 0.4/(2\pi M_{\text{final}})$, with the exact value depending on the remnant’s spin. For this to lie within the LISA band, $[10^{-4}, 10^{-1}]$ Hz, the mass of the remnant must therefore fall within the range $[10^8, 10^5]$ M$_\odot$. Ringdown modelling will be particularly important for MBHBs at the higher end of the mass spectrum, since their inspiral will take place outside the LISA band and contribute little to the SNR [1341].
4.4.3. *Status* At linear order in BHPT and for most of the duration of the signal, the ringdown is well described by a superposition of QNMs. For a given BH mass and spin, the Kerr QNFs were first computed numerically in [1334]. The most accurate and reliable numerical method is Leaver’s continued fraction approach [1342–1344] (see Mathematica notebook and data at [1345]) and in particular its spectral refinement [1156, 1346] (implemented in the Python package `qnm` [1347]). Kerr QNFs can also be estimated analytically using WKB (Wentzel–Kramers–Brillouin) or WKB-inspired approximations [1348–1353]. QNM modeling should be performed in a certain preferred BMS frame [563, 1354]. Afterwards it can be transformed to a different desired reference frame, e.g., a frame co-precessing with the binary [763, 1355], a frame set prior to a merger kick [1356], or the post-Newtonian preferred frame appropriate to modeling the early inspiral [563].

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

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

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

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

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

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

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

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

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

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

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

Coordinators: Tanja Hinderer, Geraint Pratten

Contributors: S. Akcay, A. Antonelli, S. Bernuzzi, A. Buonanno, J. Garcia-Bellido, A. Nagar, L. Pompili

#### 4.5.1. General description
The effective-one-body (EOB) approach was originally introduced in Refs. [755, 1100, 1244, 1409, 1410] with the aim of providing GW detectors with semianalytic waveform models for the entire coalescence of compact-object binaries (i.e., the inspiral, plunge, merger and ringdown), resumming PN information around the strong-field test-body limit. Since the breakthrough in NR in 2005 [514–516], the EOB framework has incorporated information from the NR simulations, thus producing highly-accurate waveform models for GW observations (e.g., see the review articles [800, 1411] and discussion below). Over the years, the EOB framework has been extended to the scattering problem, and has incorporated analytical information from other methods, such as the PM 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_1m_2/M$ is the reduced mass and $M = m_1 + m_2$ is the total mass. More specifically, the real or EOB Hamiltonian, $H_{\text{EOB}}$, is related to the effective Hamiltonian, $H_{\text{eff}}$, by [755]

$$H_{\text{EOB}} = M\sqrt{1 + 2\nu\left(\frac{H_{\text{eff}}}{\mu} - 1\right)}. \qquad (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, 916, 970, 983]. 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{dr}{dt} = \frac{\partial H_{\text{EOB}}}{\partial \boldsymbol{p}}, \qquad \frac{dp}{dt} = -\frac{\partial H_{\text{EOB}}}{\partial \boldsymbol{r}} + \boldsymbol{\mathcal{F}}, \qquad \frac{d\boldsymbol{S}_{1,2}}{dt} = \frac{\partial H_{\text{EOB}}}{\partial \boldsymbol{S}_{1,2}} \times \boldsymbol{S}_{1,2}, \qquad (16)$$

where $r$ and $p$ are the canonical variables, notably the relative position and momentum,
<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

respectively, $\mathcal{F}$ is the RR force, and $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]}, \quad (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–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].

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. [1275, 1424–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_{\ell m}^{\text{NQC}}, \quad (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 [1275, 1426, 1427]. 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_{\ell m}^{\text{NQC}}$ 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}^{\text{merger-RD}}(t)$, to the inspiral-plunge one, $h_{\ell m}^{\text{insp-plunge}}(t)$, at a suitable matching time $t = t_{\text{match}}$, around the peak of the EOB orbital frequency (which approximately corresponds to merger time), that is [1100]:

$$h_{\ell m}(t) = h_{\ell m}^{\text{insp-plunge-merger}}(t) \Theta(t_{\text{match}}^{\ell m} - t) + h_{\ell m}^{\text{ringdown}}(t) \Theta(t - t_{\text{match}}^{\ell m}), \quad (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 [575, 1100, 1126, 1247, 1251, 1410, 1424, 1425, 1434–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 [599, 1440].

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, 756, 1428, 1441]) and **TEOBResumS** (e.g., see Refs. [760, 1415, 1442, 1443]),
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 [917, 1005, 1444–1447], periastron advance [1114, 1115, 1448], scattering
angles [720, 986, 987, 1449] 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 [1420, 1450–1453], by
building reduced-order models (ROMs) [1429, 1454–1463], implementing the quasi-
adiabatic approximation [760, 1420, 1452, 1464], and making use of machine-learning
methods [1461, 1463, 1465, 1466]. 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 [575, 805, 1100, 1409, 1423, 1434-1436, 1468-1470] and spinning [754, 756, 760, 1244, 1251, 1410, 1413-1416, 1418, 1419, 1421, 1424, 1425, 1428, 1429, 1437, 1441, 1442, 1471-1476] binaries. Furthermore, orbital eccentricity [714, 752, 757-759, 1477, 1478] and matter effects [1446, 1479-1484] (see Sec. 6.5), as well as information from PM [915-917, 985-987] and conservative GSF information [505, 984, 1006, 1009, 1430, 1485-1489] 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 the two main EOB families in some more detail.

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
<table><caption>**Table 10.** Progress in the development of SEOBNR 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 SEOBNRv4PHM_surr and 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).</caption><thead><tr><th>family</th><th>waveform model</th><th>spins</th><th>ecc./hyper.</th><th>NR &amp; Teukolsky calib. region</th><th>CP-frame modes</th></tr></thead><tbody><tr><td rowspan="3">1<sup>st</sup></td><td>*EOBNRv1 [1434]</td><td></td><td></td><td>q &le; 4</td><td></td></tr><tr><td>EOBNRv2 [1469]</td><td></td><td></td><td>q &le; 3</td><td>(2, 2)</td></tr><tr><td>SEOBNRv0 [1472]</td><td></td><td></td><td>q = 1, &chi;<sub>1,2</sub> = &plusmn;0.4</td><td></td></tr><tr><td rowspan="5">2<sup>nd</sup></td><td>*EOBNRv2HM [1435]</td><td></td><td></td><td>q &le; 6</td><td>(2, 2), (2, 1), (3, 3), (4, 4), (5, 5)</td></tr><tr><td>*EOBNRv2HM_ROM [504]</td><td></td><td></td><td>q &le; 6, &chi;<sub>1,2</sub> = &plusmn;0.4</td><td></td></tr><tr><td>*SEOBNRv1 [1424]</td><td></td><td></td><td>NR: q &le; 8, |&chi;<sub>1,2</sub>| &le; 0.98</td><td>(2, 2)</td></tr><tr><td>*SEOBNRv1_ROM [1455]</td><td>✓</td><td></td><td>Teuk.: q = 1000, |&chi;<sub>1,2</sub>| &le; 0.99</td><td>(2, 2), (2, 1)</td></tr><tr><td>*SEOBNRv2 [1425]</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">3<sup>rd</sup></td><td>*SEOBNRv2_ROM [1456]</td><td></td><td></td><td></td><td></td></tr><tr><td>SEOBNRv2P [1251]</td><td>✓✓</td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv3P [1425, 1437]</td><td>✓✓</td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv4 [754]</td><td>✓</td><td></td><td></td><td>(2, 2)</td></tr><tr><td rowspan="7">4<sup>th</sup></td><td>*SEOBNRv4_ROM [754]</td><td></td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv4P [756]</td><td>✓✓</td><td></td><td></td><td>(2, 2), (2, 1)</td></tr><tr><td>*SEOBNRv4HM [1428]</td><td></td><td></td><td>NR: q &le; 8, |&chi;<sub>1,2</sub>| &le; 0.995</td><td></td></tr><tr><td>*SEOBNRv4HM_ROM [1459]</td><td>✓</td><td></td><td>Teuk.: q = 1000, |&chi;<sub>1,2</sub>| &le; 0.99</td><td>(2, 2), (2, 1), (3, 3), (4, 4), (5, 5)</td></tr><tr><td>*SEOBNRv4PHM [756]</td><td></td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv4PHM_surr [1460]</td><td>✓✓</td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv4EHM [758]</td><td>✓</td><td>✓</td><td></td><td>(2, 2), (2, 1), (3, 3), (4, 4)</td></tr><tr><td>*SEOBNRv4PHM_4dq2 [1463]</td><td>✓✓</td><td></td><td></td><td></td></tr><tr><td rowspan="4">5<sup>th</sup></td><td>*SEOBNRv5 [1429]</td><td></td><td></td><td>NR: q &le; 20, |&chi;<sub>1,2</sub>| &le; 0.998</td><td>(2, 2)</td></tr><tr><td>*SEOBNRv5_ROM [1429]</td><td>✓</td><td></td><td>Teuk.: q = 1000, |&chi;<sub>1,2</sub>| &le; 0.99</td><td>(2, 2), (2, 1), (3, 2), (3, 3), (4, 3), (4, 4), (5, 5)</td></tr><tr><td>*SEOBNRv5HM [1429]</td><td></td><td></td><td></td><td></td></tr><tr><td>*SEOBNRv5PHM [1441]</td><td>✓✓</td><td></td><td></td><td></td></tr></tbody><tfoot><tr><td colspan="6">✓ aligned spin, ✓✓ arbitrary spin orientations.<br>* implemented in LALsuite, and for the latest version (v5) in the open-source Python package pySEOBNR [1490].</td></tr></tfoot></table>

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 [984, 1009, 1430, 1486–1488], 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 § 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 10. 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, 1244, 1409], 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–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 [1158, 1246, 1423], 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–1419]. It included all PN § The generic name SEOBNRvnnEPHM 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).
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, 752, 1498] 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 [113, 115, 210, 1499, 1500].

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-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. [760, 1502, 1503], 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 ~ 8-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 [765, 1504, 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 [754, 1429, 1455, 1456, 1460]. 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. [1009, 1127, 1425, 1486, 1487], and to non-vacuum binary systems and gravity theories beyond GR are discussed in Sec. 6.5.

4.5.4. The TEOBResumS waveform models 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 [1247, 1423], the systematic resummation of EOB potentials via Padé approximants [1247, 1413, 1414, 1436, 1507], NQC corrections [1245, 1246, 1436, 1471, 1473], 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–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 [1128, 1245, 1246, 1249, 1250] 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,
<table><caption>**Table 11.** 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.</caption><thead><tr><th>Physical content</th><th>GIOTTO</th><th>DALI</th></tr></thead><tbody><tr><td>Analytic information</td><td>A(r)<br>D(r)<br>Q(r, p<sub>r*</sub>)<br>G<sub>S</sub>, G<sub>S*</sub><br>r<sub>c</sub></td><td>Pseudo 5PN, resummed<br>3PN, resummed<br>3PN<br>3.5PN<br>NLO</td><td>Pseudo 5PN, resummed<br>5PN, resummed<br>5PN (local)<br>3.5PN<br>NLO</td></tr><tr><td>NR information</td><td>a<sub>c</sub><sup>6</sup><br>c<sub>N3LO</sub><br>NQCs<br>Ringdown<br>BBH NR Validation region</td><td>Effective parameter in A(r)<br>Effective parameter in G<sub>S</sub>, G<sub>S*</sub><br>Ensure correct transition between plunge and merger<br>Phenomenological model, quasi-circular<br>q ≥ 10 and test-mass with |χ<sub>1,2</sub>| ≥ 0.99; 10 &lt; q &lt; 128 no-spins</td><td></td></tr><tr><td>Spins</td><td>Aligned<br>Precessing</td><td>✓<br>✓</td><td>✓<br>-</td></tr><tr><td>Orbital dynamics</td><td></td><td>Circular</td><td>Generic (bound & open)</td></tr><tr><td>Subdominant CP modes</td><td>Inspiral to merger<br>Merger/ringdown</td><td>(ℓ, |m|) ≤ 8<br>(ℓ, |m|) = (2, 1), (2, 2), (3, 2), (3, 3), (4, 2), (4, 4), (5, 5)</td><td></td></tr></tbody></table>

BHNS)[1484, 1508, 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 [720, 1512, 1513]. The key features of the models are summarized in Table 11 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 [505, 1489], and work on scalar-tensor gravity (building upon Refs. [1514, 1515]) is also currently in progress [1516–1518].

The structure of the spin-aligned effective TEOBResumS Hamiltonian is

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

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

$$H_{\text{eff}}^{\text{SO}} = G_{S}S + G_{S*}S_{*}, \qquad (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 [1413, 1414, 1519]. 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^3LO$ 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” [760, 1476, 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. [1427, 1443]. 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–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, 1461, 1475, 1508] 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.
<figure><img src="image_11.png" /><figcaption>**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 $l = 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].</figcaption></figure>

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:

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

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

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 [522, 728, 753, 761, 762, 1525–1537], 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 $l = |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 ~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 ≈ $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, 761, 1502, 1529–1532], 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, 761, 1529-1532] 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:

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

(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.

(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:

* PhenomPv2 [1529] updates the PhenomP model for precession.

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

* 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].

* 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 [761, 762, 1536, 1537] 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:

* The number of NR waveforms used for calibration is increased from 20 to ~400.

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

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

* Both spin degrees of freedom are calibrated to numerical data.

* 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 [761, 1537]. 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
<table><caption>**Table 12.** Progress in the development of phenomenological waveform models in frequency domain (FD) and time domain (TD).</caption><thead><tr><th>Waveform Family</th><th>Domain</th><th>Waveform Model</th><th>Spins</th><th>Mode Content</th><th>Eccentricity</th><th>Calibration Region</th></tr></thead><tbody><tr><td rowspan="3">1st generation</td><td rowspan="3">2nd generation</td><td>IMRPhenomA</td><td>✗</td><td rowspan="3">(2,±2)</td><td rowspan="10">no</td><td>q ≤ 4</td></tr><tr><td>IMRPhenomB</td><td>✓</td><td rowspan="2">NR calibration:<br>q ≤ 4, |χ<sub>1/2</sub>| ≤ 0.75<br>|χ<sub>1/2</sub>| ≤ 0.85 (for q = 1)</td></tr><tr><td>IMRPhenomC</td><td>✓</td></tr><tr><td rowspan="7">3rd generation</td><td rowspan="7">FD</td><td>IMRPhenomP</td><td>✓✓</td><td rowspan="2">CP</td><td rowspan="2">NR calibration:<br>q ≤ 18, |χ<sub>1/2</sub>| ≤ 0.85<br>-0.95 ≤ χ<sub>1/2</sub> ≤ 0.98<br>(for q = 1)</td></tr><tr><td>IMRPhenomD</td><td>✓</td></tr><tr><td>IMRPhenomPv2</td><td>✓✓</td><td rowspan="2">CP</td></tr><tr><td>IMRPhenomPv3</td><td>✓✓</td></tr><tr><td>IMRPhenomHM</td><td>✓</td><td rowspan="2">(2,±2), (2,±1), (3,±3),<br>(4,±3), (4,±4)</td><td rowspan="4">in development</td><td rowspan="4">NR calibration:<br>q ≤ 18, |χ<sub>1/2</sub>| ≤ 0.99<br>Teukolsky calibration:<br>q ≤ 1000</td></tr><tr><td>IMRPhenomPv3HM</td><td>✓✓</td></tr><tr><td>IMRPhenomXAS</td><td>✓</td><td rowspan="2">(2,±2)</td></tr><tr><td>IMRPhenomXP</td><td>✓✓</td></tr><tr><td rowspan="6">4th generation</td><td rowspan="6">TD</td><td>IMRPhenomXHM</td><td>✓</td><td rowspan="2">(2,±2), (2,±1), (3,±2),<br>(3,±3), (4,±4)</td><td rowspan="6">in development</td><td rowspan="6">NR calibration:<br>q ≤ 18, |χ<sub>1/2</sub>| ≤ 0.99<br>Teukolsky calibration:<br>q ≤ 1000</td></tr><tr><td>IMRPhenomXPHM</td><td>✓✓</td></tr><tr><td>IMRPhenomT</td><td>✓</td><td rowspan="2">(2,±2)</td></tr><tr><td>IMRPhenomTP</td><td>✓✓</td></tr><tr><td>IMRPhenomTHM</td><td>✓</td><td rowspan="2">(2,±2), (2,±1), (3,±3),<br>(4,±4), (5,±5)</td><td rowspan="2">in development</td><td rowspan="2">NR calibration:<br>q ≤ 18, |χ<sub>1/2</sub>| ≤ 0.99<br>Teukolsky calibration:<br>q ≤ 1000</td></tr><tr><td>IMRPhenomTPHM</td><td>✓✓</td></tr></tbody><tfoot><tr><td colspan="7">✗ no spins ✓ spins aligned with orbital angular momentum ✓✓ precessing spins CP mode content in co-precessing frame</td></tr></tfoot></table>

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. 12.

### 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:

(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.

(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:

(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.)

(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.

(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.

(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 [1540, 1544, 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.

### 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 $(l, 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 [589, 590, 1548]. It is often advantageous to further decompose the data into symmetric and asymmetric mode combinations [589, 590, 1548]. 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.
<figure><img src="image_12.png" /><figcaption>**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 $\lambda_0$, that is not in the training set. To compute the waveform for this value, each fit is evaluated at $\lambda_0$ (the yellow diamonds), and then the reduced-order representation is used to reconstruct the decompressed waveform (the dotted black line).</figcaption></figure>

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 [1142, 1156]. 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–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
[1454, 1455, 1559, 1560]. 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 [1462, 1463, 1562, 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 [1556–1560, 1565]. Since then, ROM techniques have been further developed and refined for various LVK (hence comparable-mass) applications: EOB [1454–1456, 1459, 1460], NR [589, 590, 764, 1504, 1548], multiple subdominant modes [589, 590, 764, 1454, 1459, 1504, 1548], frequency-domain waveforms [1455, 1456, 1459], time-domain waveforms [589, 590, 764, 1454, 1504, 1548], generically precessing systems [589, 590, 764, 1548], neutron star inspirals with tidal effects [1457, 1458, 1566], and
waveforms with ~ $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 [114, 1500, 1567-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 [215, 487, 488]. 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 [766, 1117] 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,000M 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 ~ $10^5$ constituent modes. Fast frameworks for EMRI trajectory generation already exist [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 (~ $\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 (~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] (~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-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 (~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 (~2000x) [1110, 1160]. 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 ~ $10^6 - 10^7$ time points, each with a combination of ~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 ~ 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. [275, 1575–1582]. 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. [275, 1578]. 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. [614, 630, 638, 1583–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 [614, 638, 1583, 1594]. 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 [1591, 1592, 1598–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* [630, 1584, 1586, 1587, 1615–1617], in which the spacetime metric and extra fields are expanded around GR, order-by-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 [1593, 1601, 1618–1620]. 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 [630, 1584, 1586, 1587, 1593, 1612, 1615–1617, 1619, 1621, 1622], 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–1625], Einstein dilaton Gauss-Bonnet gravity (EdGB) [334, 403, 1626, 1627] or more generic Einstein scalar Gauss-Bonnet theories (ESGB) [1628–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–1631] (see also [1622, 1632–1637] 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 [630, 1621]. 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 full BBH parameter space, build a surrogate model for parameter estimation [590], and perform model-dependent tests of GR with gravitational

*Sometimes 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.*
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, 1620, 1640, 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 [1585, 1647-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, 355-357, 364, 366-368, 373, 1654-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 [374, 375, 378, 379, 1661-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, 438, 1664-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-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 [375, 1663, 1683, 1690, 1691]. 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–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 [1697, 1704–1706]. 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–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 M f)^{a_i}]e^{i\beta_j(\pi M f)^{b_j}}, \quad (22) $$

where $\tilde{h}_{GR}(f)$ is the frequency domain waveform of GR, $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. [1575, 1719–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 [1718, 1729], 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 [64, 1735, 1736]. Such an environment can affect the binary dynamics through dynamical friction [1737], accretion and the gravitational pull of the environment itself [64, 1735, 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 [372, 374, 380–383, 1740–1743] showing that a plethora of signatures, including tidally induced resonances [380, 382], floating and kicking orbits [382, 1740, 1744, 1745] and ionization of the cloud [1742, 1746–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 [377, 378, 1749], 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–470], which at leading order enters as a 2PN correction to the GW phase due to the object’s quadrupole moment [469, 1750]. 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 [830, 1156, 1751–1753]. Absorption at the horizon introduces a 2.5PN (4PN) × log*ν* 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–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 [612, 1335, 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 [606, 612, 1335, 1336, 1358, 1364, 1367–1369, 1371, 1374, 1762–1770]. 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 [1354, 1361, 1363, 1368–1370, 1373, 1376, 1393, 1438, 1771–1773].

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–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–1783], or they are given in the form of numerical solutions of the field equations (e.g. [359, 360, 1784–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. [64, 459, 1782, 1790–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. [1777, 1801–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 [1354, 1361–1363, 1369, 1370, 1373, 1375, 1376, 1393, 1771–1773], 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–1811], wormholes [437, 1812, 1813], quantum structures [438, 1814–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 [459, 1820, 1821]. 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 [436, 1329, 1797, 1817]. The formation of an exotic ultracompact object may
also lead to the emission of GW echoes in the post-merger signal [434, 437] (see also
Sec. 2.7). Several phenomenological GW templates for echoes have been developed based
on standard GR ringdown templates with additional parameters [433, 1822, 1823], the
superposition of sine-Gaussians [1824], hybrid methods that put together information
from perturbation theory and the pre-merger orbital dynamics [1825–1827], using the
close-limit approximation [1828], and analytical models that explicitly depend on the
physical parameters of the ECO [456, 1829–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, 1744, 1792, 1832–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–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 [1069, 1098, 1240, 1241, 1301, 1844–1852].

**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 [64, 376, 1736, 1853, 1854] (accretion and the usual gravitational pull of the dark matter can be significant, but are in general subleading effects [64, 1735, 1736, 1855]). 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–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, 1289, 1859, 1862–1864]. 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 [372, 380, 383–385, 1735, 1740, 1745–1748, 1853, 1855, 1859, 1865–1867]. 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. [1290, 1291, 1868] 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. [1290, 1291, 1868, 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 [332, 335, 342, 343, 365, 369–371, 387, 389, 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 [343, 365, 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}$ [99, 105], 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 [104, 468, 1871–1874]. 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 [471, 472, 1282, 1758, 1877–1883]. Because bumpy BHs break the symmetries of Kerr that yield integrable motion [1884], orbits in these spacetimes are generically chaotic [1885–1888] and subject to prolonged resonances [1318, 1887–1889]. Most of the work considering EMRIs around spacetimes with a multipolar structure different from Kerr, constructed approximate “kludge” waveforms [468, 488, 1874, 1890, 1891] generated considering geodesics in a perturbed Kerr spacetime with parameters evolved using post-Newtonian equations [105, 107, 1233, 1873, 1885]. 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 [457, 1756, 1893–1896]. In addition, the modes of the ECO could also be excited during the inspiral of the small compact object [457, 1807, 1894, 1895, 1897–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 [457, 1894, 1895]. 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 [604, 605, 1718, 1902, 1903], 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 [604–606, 1764, 1904, 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–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. [1038, 1442, 1446, 1452, 1475, 1481, 1484, 1909–1911], with the state-of-the-art models incorporating the known adiabatic gravitoelectric and -magnetic tidal terms up to $\ell = 8$ in different resummation schemes [1481, 1910–1912], dynamical tides from the fundamental modes [1479, 1480, 1913] 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 [1479, 1480], 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 [1533, 1917],
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 [1918–1920], and for the tidal disruption in BH-neutron-
star binaries [1483, 1484, 1921]. 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 [604, 605, 1902,
1924, 1925]. 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 [282, 283, 1928–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 [263, 1930] 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, 302, 305, 1936].

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–1941]. Gravitational backreaction probably smoothens the string at small scales [1942–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 closely approach, but do not strictly reach, the cusp configuration. These pseudocusp burst waveforms are of lower amplitude than true cusps, and follow the same power spectrum decay of $f^{-4/3}$ at high frequencies. Accounting for them in a search can lead to an enhancement in the number of expected burst events.
**6.7.4. Beyond Nambu–Goto for cusps** Gauge theory strings with energy scale $\mu_S^{1/2}$ have a characteristic core width $\sim \mu_S^{-1/2}$ whereas classical superstrings are always well-represented by the zero thickness Nambu–Goto limit. When a cusp forms the string instantaneously doubles back upon itself. Ref. [1948] predicts that a region of string around the cusp, of the order of the string thickness, will annihilate into particles due to overlap. For strings with large tensions (near that of the GUT scale) this modification to the cusp shape is at a small enough scale that it does not meaningfully affect the cusp waveform. However, as tension decreases the physical scale increases. It may prove useful to compare the waveform of the gauge theory cusp to the superstring cusp. This can be done by parametrically varying the core width in Abelian-Higgs simulations.

A final assumption about the standard cusp waveform model is that it relies only on properties universal to all strings. However, superstrings possess additional features that could modify their cusp waveforms. Ref. [298] suggested modifications to the amplitude based on the dimensionless loop-length parameter $\alpha$. For example, cusps very near to Y-junctions will have their burst amplitude modified by a correction factor that depends on the nearby string geometry [271], while the extra degrees of freedom on superstrings (extra dimensions) can lead to a strong suppression of both cusp rate and burst amplitude [1949].

**6.7.5. Cosmic strings and GW production in numerical relativity** While in most cases, perturbative methods (e.g., Nambu–Goto strings coupled to perturbative gravity [1943]) are sufficient, situations may arise when perturbative methods are insufficient and gravitational backreaction becomes of $O(1)$, even if the string tension $G\mu_S$ is small. For example, a string with a sharp kink can easily lead to a locally strong field situation [1950]. In such strong field limits, the full machinery of NR must be employed. The construction of full NR solutions for cosmic strings was recently achieved by [1951, 1952], where stable dynamical solutions for an Abelian-Higgs string minimally coupled to gravity was found. This opens the door to not just the detailed exploration of the string dynamics, but also the characterization of the GW signals that can then be searched for.

The Abelian-Higgs model minimally coupled to gravity can be recast into the standard BSSN/CCZ4 formalism, and solved using finite differences (we refer the reader to [1951] for details). Presently, solutions for infinite static strings, planar loops [1951] and single traveling wave solutions [1953] have been constructed. Once a string configuration is constructed and successfully evolved, the work required to extract useful GW signals from them follows very similarly to that of any standard NR simulations, see Sec. 4.1.

As an example, in [1952], the gravitational waveform emitted during the BH formation of a collapsing circular cosmic string was computed, and it was shown that the waveform is dominated by the BH formation and ringdown phase, with the primary contribution being the $l = 2, m = 0$ mode (as opposed to quasi-circular compact binary merger
signals which are dominated by the $ℓ = 2, m = 2$ mode). Intriguingly, due to the large
asymmetric ejection of material in the formation process, a large GW memory was also
seen.

The future prospects for using NR to compute GW signals from such strong gravity
events is very promising. For example, a full general relativistic treatment of the
GW emission from cosmic string cusps and kinks is long overdue, but is now within
reach. Presently, the main obstacle to explore more complicated string configurations
is the need for the construction of accurate initial conditions – a usual problem in any
NR endeavor, here made more difficult by the presence of fundamental matter fields.
Nevertheless, this area is ripe for further exploration.
7. Afterword

LISA will usher in an era of millihertz GW astronomy that will open a window to new
classes of sources and offer unprecedented opportunities to probe our understanding of
the universe. Integral to this scientific vision are predictions of the source waveforms,
which are necessary for the detection and interpretation of GW events.

This white paper has analyzed the question of waveform preparedness in the LISA
era. It has reviewed the modelling requirements from astrophysical and data analysis
perspectives, and has analyzed the status of the main approaches to modelling waveforms
from compact binaries both in vacuum GR and in the presence of environmental
influences or effects from new physics beyond GR and the Standard Model. This white
paper also provides guidance on new developments that are needed for the different
approaches, both in terms of modelling accuracy and the covered parameter space, so
that gravitational wave models will be ready for the LISA era.

The experience with ground-based GW detectors has to some extent prepared the
community. For example, several approaches that predict the waveforms for comparable
mass black hole binaries can be applied directly to massive black hole binaries due to the
simple mass scaling of GR. However, these models are not currently sufficient to fully
realize LISA's science goals. In particular they insufficiently cover the parameters of
eccentricity, highly asymmetric mass-ratios, and high spin magnitudes. In addition, the
possibility of very strong massive black hole binary signals places additional demands
on the accuracy of these models, well beyond what can currently be achieved.

For extreme-mass-ratio inspirals, leading-order (adiabatic) models that are sufficient to
enable some LISA science are becoming available, and these models can be evaluted
fast enough for use in LISA data analysis. In order to enable the full spectrum of
LISA science with EMRs post-adiabatic models are needed. These have recently been
developed for the simplest orbital configurations (non-spinning, quasi-circular). It also
looks likely that these models can cover much of the intermediate mass ratio parameter
space. Significant work remains to extend these post-adiabatic models to cover the full
parameter space of precessing and eccentric parameter binaries where both components
are spinning.

The modelling of sources in theories beyond GR is still in its infancy with significant
effort needed. For example, the majority of computations have focused on higher
derivative gravity, that are most interesting in the high curvature regime, or on
Horndeski models. Moreover, simulations of compact binaries in extensions of GR are,
at the time of writing this white paper, proof-of-concept. That is, they have been done
for nearly equal-mass binaries, typically of non spinning black holes, and at relatively
low resolution. This is comparable to the state of numerical relativity in GR about 15
years ago (i.e., few years after the breakthrough in 2005). Likewise, perturbative or
PN calculations have only started to go beyond the leading order in the correction. To
perform theory-specific tests of gravity, more work is needed, both in terms of modelling
accuracy and in covering (parts of) the parameter space of black hole masses, spins and
additional theory-specific parameters.

The good news is that in all approaches there are clear ideas of what needs to be done
in order to reach the modelling requirements for LISA. Doing so will require a concerted
effort over the next decade. Given sufficient resources, we are positive that all goals can
be achieved.

We particularly recommend that the community conduct more investigations into the
waveform standards and waveform interface. While lessons have been learned from
ground-based detectors, LISA will offer unique challenges that will need to be addressed
by the waveform and data analysis communities working together.
Appendix A: Descriptions of NR codes

Description of some numerical relativity codes that focus, in particular, on modelling
LISA sources. A more extensive list is given in Table 8.

**BAM** BAM is a modular code framework initially developed at the University of Jena and is now maintained and further developed by the CoRe (Computational Relativity) Collaboration [582, 616–618]. BAM employs finite differences for the discretization of spacetime fields, and high-resolution shock-capturing methods for the evolution of general relativistic hydrodynamics variables. It features automated code generation with Mathematica, adaptive mesh refinement (AMR) in space and time, and parallelization with MPI and OpenMP for large scale HPC, recently at 70% on up to 50k cores. BAM is employed for the simulation of compact binary mergers, e.g., [524, 538, 1910, 1954]. A major focus is on binary black hole, binary neutron star, and black hole - neutron star systems [1955], plus more exotic compact objects or merger scenarios [1685]. Recent updates allow the study of radiation hydrodynamics [1956], magnetohydrodynamics, and the usage of an entropy-viscosity limiter for the flux computation [1957]. BAM has contributed numerical data for the development and validation of e.g. Phenom-type and EOB waveform models. Several hundred binary neutron star merger simulations are available as part of the existing CoRe database [619, 620].

**BAMPS** The Jena-Lisbon collaboration is developing the code bamps, a new and highly scalable NR code for exascale applications using pseudo-spectral and DG methods [621, 622]. A key feature is hp-refinement for spectral element methods [623], which allows highly scalable and accurate simulations, currently with a focus on smooth fields. To date, bamps has already been employed for the study of GW collapse [622, 1958], simulations with scalar fields for collapse and for boson stars are in development. Further applications include simple neutron star spacetimes [621], and studies of new techniques regarding the dual foliation framework [1959].

**Dendro-GR** Dendro-GR [626, 1960, 1961] uses an unstructured grid with a localized, wavelet-based refinement scheme to evolve the BSSN formulation of the Einstein equations using moving punctures. In test simulations of a BBH with a mass ratio of 10 : 1, Dendro-GR shows good scaling to over $10^5$ cores.

**Einstein Toolkit** The Einstein Toolkit [628, 629] is an open-source cyberinfrastructure for computational astrophysics, with more than 300 registered users worldwide. It implements Einstein's equations, or extensions thereof, in the BSSN or CCZ4 formulation together with the moving puncture gauge. It employs the method of lines, in which spatial derivatives are implemented as up to 8th order finite differences together with a collection of direct time integration techniques (e.g., fourth order Runge-Kutta integrator). The Einstein Toolkit is based on the Cactus computational toolkit [1962] and the Carpet boxes-in-boxes AMR package [1963,
[1964]. It also provides a multipatch infrastructure consisting of Cartesian boxes-in-boxes meshes in the region where the binaries evolve (typical size ~50...100M), and a spherical outer region in the wave zone (typical radius ~1000...2000M) provided by the Llama code. Kreiss-Oliger dissipation is employed to reduce high frequency noise at refinement boundaries. The toolkit uses hybrid MPI/OpenMP parallelization. The Einstein Toolkit has been a broad community project with a multitude of software that is based on its infrastructure or that has been developed within its framework. The additional software is indicated by asterisks in Table 8.

**Einstein Toolkit – CarpetX** The Einstein Toolkit consortium is developing a new AMR driver, CarpetX [1965], for the Einstein Toolkit. CarpetX leverages the AMReX framework [1966] to provide scalable adaptive mesh refinement for physics modules. CarpetX will provide: efficient support for both CPUs and GPUs; parallelization via MPI and OpenMP; SIMD vectorization; scalable I/O based on the ADIOS2 file format [1967] and the openPMD metadata standard; and more. CarpetX offers block-structured AMR based on local error estimates, exact conservation across mesh interfaces, higher order prolongation operators, and scalable elliptic solvers based on PETSc [1968]. CarpetX is available as open source and was used in production in Ref. [1969].

**GR-Athena++** Another example is GR-Athena++ [652] based on the astrophysical (radiation) magnetohydrodynamics code Athena++ [1970]. GR-Athena++ leverages on Athena++'s oct-tree AMR and exploits hybrid parallelism at different levels. The standard distributed approach is augmented by a dynamical tasklist for the procedures in each mesh block that overlaps calculation and communication (thus mitigating some of the overhead associated with high order finite differencing). In-core vectorization and other optimization techniques of the basic kernels [1971] are also employed. This approach leads to a high parallel efficiency, with strong scaling efficiencies above 95% for up to $10^4$ CPUs and excellent weak scaling is up to $10^5$ CPUs measured in a production BBH setup with AMR.

**GRChombo** GRChombo ([653–655]) is an open-source NR code built using an optimized version of the publicly available library Chombo [1972] developed at LBNL. It is based on established methods of solving the Einstein equations, but the highly flexible AMR capability and templating over physics classes makes it suited to strong gravity problems in beyond-GR scenarios (e.g. Horndeski gravity [1973], cosmic strings [1952]), which have been key targets of the code. GRChombo is written entirely in C++14, using hybrid MPI/OpenMP parallelism and vector intrinsics (in particular, the evolution calculations are explicitly vectorised) to achieve improved performance on the latest architectures. It has good strong scaling up to around 4,000 cores for typical binary BH problems, at which point it becomes limited by problem size as with most traditional codes.

The next release of Chombo is designed to be performance portable to different heterogeneous architectures (including GPUs), which GRChombo should be able
to leverage to improve scaling further.

**Illinois GRMHD** The Illinois GRMHD evolution code [659, 660, 1974] solves the Einstein field equations via the BSSN scheme with puncture gauge conditions, utilizes a conservative high-resolution, shock-capturing (HRSC) scheme for the matter and adapts the Carpet infrastructure to implement AMR. The code solves the magnetic induction equation by introducing a magnetic vector potential, which guarantees that the B-field remains divergence-free for any interpolation scheme used on refinement level boundaries. This formulation reduces to a standard constrained transport (CT) scheme on uniform grids. An improved Lorenz gauge condition for evolving the vector potential in AMR without the appearance of spurious B-fields on refinement level boundaries is implemented. This gauge exhibits no zero-speed modes, enabling spurious magnetic effects to propagate off the grid quickly. The GRMHD code recently added a radiation module built to handle transport (photons or neutrinos) via the M1 moment formalism. A version of the Illinois GRMHD code has been released as an open-source module that has been ported to the Einstein Toolkit, where it has been documented and rewritten to make it more user-friendly, modular and efficient. The code has been used to simulate numerous scenarios related to LISA sources (e.g., BHBH mergers, BHBH mergers in gaseous clouds and magnetized disks, etc.) and their associated gravitational waveforms.

**MHDuet** The publicly available code MHDuet [663, 664, 1975, 1976] is generated by the open-source platform SIMFLOWNY to run under the SAMRAI infrastructure, which provides the parallelization and the adaptive mesh refinement. There are different versions of the code, either to study alternative gravity theories or to include large-eddy-simulations and neutrino transport.

**SACRA** The SimulAtor for Compact objects in Relativistic Astrophysics (SACRA) code solves Einstein's equation with the BSSN-puncture formulation together with a Z4c constraint propagation prescription and the relativistic hydrodynamics equation with the HRSC scheme [666, 667]. It implements the box-in-box conservative AMR with parallelization by MPI and OpenMP. The code has been used to simulate compact binary mergers such as NSNS mergers, BHNS mergers, and BHBH mergers. The NSNS waveform catalog SACRA Gravitational Wave Data Bank is available at [671]. The latest version of SACRA has a couple of branches, such as SACRA-TD, to simulate the tidal disruption of an ordinary star by the SMBH [1977], or SACRA-MG to simulate a compact object in scalar-Gauss-Bonnet theory [668] and scalar-tensor theory [669]. SACRA also has a version to solve high-dimensional Einstein's equation [670]. In addition, it features modules for neutrino radiation transport [1978], MHD [1979, 1980] and viscous hydrodynamics [1981].

**SACRA-SFS2D** is a viscous-radiation hydrodynamics [672] and radiation magnetohydrodynamics code [673] in full general relativity. These codes can be applied to the