# Waveform Modelling for the Laser Interferometer Space Antenna

Niaesh Afshordi, Sarp Akçay, Pau Amaro Seoane, Andrea Antonelli, Josu C Aurrekoetxea, Leor Barack, Enrico Barausse, Robert Benkel, Laura Bernard, Sebastiano Bernuzzi, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Niaesh Afshordi, Sarp Akçay, Pau Amaro Seoane, Andrea Antonelli, Josu C Aurrekoetxea, et al. "Waveform Modelling for the Laser Interferometer Space Antenna." Living Reviews in Relativity, 2025, 28 (9), 10.1007/s41114-025-00056-1. [hal-04285017]

HAL Id: hal-04285017

https://hal.science/hal-04285017v1

Submitted on 4 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Waveform Modelling for the Laser Interferometer Space
Antenna

LISA Consortium Waveform Working Group

**Abstract.** LISA, the Laser Interferometer Space Antenna, will usher in a new era in gravitational-wave astronomy. As the first anticipated space-based gravitational-wave detector, it will expand our view to the millihertz gravitational-wave sky, where a spectacular variety of interesting new sources abound: from millions of ultra-compact binaries in our Galaxy, to mergers of massive black holes at cosmological distances; from the beginnings of inspirals that will venture into the ground-based detectors' view to the death spiral of compact objects into massive black holes, and many sources in between. Central to realising LISA's discovery potential are waveform models, the theoretical and phenomenological predictions of the pattern of gravitational waves that these sources emit. This white paper is presented on behalf of the Waveform Working Group for the LISA Consortium. It provides a review of the current state of waveform models for LISA sources, and describes the significant challenges that must yet be overcome.
Niayesh Afshordi,¹ Sarp Akçay,² Pau Amaro Seoane,³,⁴,⁵,⁶ Andrea Antonelli,⁷
Josu C. Aurrekoetxea,⁸ Leor Barack,⁹ Enrico Barausse,¹⁰ Robert Benkel,¹¹
Laura Bernard,¹² Sebastiano Bernuzzi,¹³ Emanuele Berti,⁷ Matteo Bonetti,¹⁴,¹⁵
Béatrice Bonga,¹⁶ Gabriele Bozzola,¹⁷ Richard Brito,¹⁸ Alessandra Buonanno,¹¹
Alejandro Cárdenas-Avendaño,¹⁹ Marc Casals,²⁰,²¹,² David F. Chernoff,²²
Alvin J. K. Chua,²³,²⁴ Katy Clough,²⁵ Marta Colleoni,²⁶ Geoffrey Compère,²⁷
Mekhi Dhesi,⁹ Adrien Druart,²⁸ Leanne Durkan,²⁹ Guillaume Faye,³⁰ Deborah
Ferguson,²⁹,³¹ Scott E. Field,³² William E. Gabella,³³ Juan García-Bellido,³⁴
Miguel Gracia-Linares,²⁹ Davide Gerosa,¹⁴,¹⁵,³⁵ Stephen R. Green,³⁶ Maria
Haney,³⁷ Mark Hannam,³⁸ Anna Heffernan,³⁹,¹,⁴⁰ Tanja Hinderer,⁴¹ Thomas
Helfer,⁴² Scott A. Hughes,⁴³ Sascha Husa,⁴⁴,³⁹ Soichiro Isoyama,²³ Michael L.
Katz,⁴⁵,¹¹ Chris Kavanagh,² Gaurav Khanna,⁴⁶ Larry E. Kidder,⁴⁷ Valeriya
Korol,⁴⁸ Lorenzo Küchler,⁹,²⁸,⁴⁹ Pablo Laguna,²⁹ François Larrouturou,⁵⁰
Alexandre Le Tiec,¹² Benjamin Leather,¹¹,² Eugene A. Lim,⁵¹ Hyun Lim,⁵²
Tyson B. Littenberg,⁴⁵ Oliver Long,¹¹,⁹ Carlos O. Lousto,⁵³ Geoffrey
Lovelace,⁵⁴ Georgios Lukes-Gerakopoulos,⁵⁵ Philip Lynch,¹¹,² Rodrigo
P. Macedo,⁵⁶ Charalampos Markakis,²⁵ Elisa Maggio,¹¹ Ilya Mandel,⁵⁷,⁵⁸ Andrea
Maselli,⁵⁹,⁶⁰ Josh Mathews,²³ Pierre Mourier,³⁹ David Neilsen,⁶¹ Alessandro
Nagar,⁶²,⁶³ David A. Nichols,⁶⁴ Jan Novák,⁶⁵ Maria Okounkova,⁶⁶ Richard
O'Shaughnessy,⁵³ Naritaka Oshita,⁶⁷,⁶⁸,⁶⁹ Conor O'Toole,² Zhen Pan,¹ Paolo
Pani,⁷⁰ George Pappas,⁷¹ Vasileios Paschalidis,⁷² Harald P. Pfeiffer,¹¹ Lorenzo
Pompili,¹¹ Adam Pound,⁹ Geraint Pratten,³⁵ Hannes R. Rüter,¹¹ Milton
Ruiz,⁷³,³¹ Zeyd Sam,⁹ Laura Sberna,¹¹ Stuart L. Shapiro,³¹ Deirdre M.
Shoemaker,²⁹ Carlos F. Sopuerta,⁴⁴,⁷⁴ Andrew Spiers,⁹ Hari Sundar,⁷⁵
Nicola Tamanini,⁷⁶ Jonathan E. Thompson,⁷⁷,³⁸ Alexandre Toubiana,¹¹ Antonios
Tsokaros,³¹,⁷⁸ Samuel D. Upton,⁵⁵,⁹ Maarten van de Meent,⁵⁶,¹¹ Daniele
Vernieri,⁷⁹ Jeremy M. Wachter,⁸⁰ Niels Warburton,² Barry Wardell,²
Helvi Witek,³¹ Vojtěch Witzany,⁸¹ Huan Yang,¹,⁴⁰ and Miguel Zilhão,⁸²
(Coordinators and contributors)

Angelica Albertini,⁵⁵,⁸¹ K. G. Arun,⁸³ Miguel Bezares,⁸⁴,⁸⁵ Alexander
Bonilla,⁸⁶ Christian Chapman-Bird,⁸⁷ Bradley Cownden,² Kevin
Cunningham,² Chris Devitt,² Sam Dolan,⁸⁸ Francisco Duque,¹¹ Conor
Dyson,⁵⁶ Chris L. Fryer,⁵² Jonathan R. Gair,¹¹ Bruno Giacomazzo,¹⁴,¹⁵,⁸⁹ Priti
Gupta,⁹⁰ Wen-Biao Han,⁹¹ Roland Haas,⁹² Eric W. Hirschmann,⁶¹
E. A. Huerta,⁹³,⁹⁴,⁹⁵ Philippe Jetzer,⁹⁶ Bernard Kelly,⁹⁷,⁹⁸,⁹⁹ Mohammed
Khalil,¹ Jack Lewis,⁹ Nicole Lloyd-Ronning,¹⁰⁰ Sylvain Marsat,⁷⁶ Germano
Nardini,¹⁰¹ Jakob Neef,² Adrian Ottewill,² Christiana Pantelidou,² Gabriel
Andres Piovano,² Jaime Redondo-Yuste,⁵⁶ Laura Sagunski,¹⁰² Leo C.
Stein,¹⁰³ Viktor Skoupý,⁸¹,⁵⁵ Ulrich Sperhake,¹⁰⁴ Lorenzo Speri,¹¹ Thomas
F.M. Spieksma,⁵⁶ Chris Stevens,¹⁰⁵ David Trestini,¹⁰⁶ and Alex Vañó-Viñuales ¹⁸
(Endorsers)

¹Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada

²School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland

³Universitat Politècnica de València, València, Spain

⁴Max Planck Institute for Extraterrestrial Physics, Garching, Germany

⁵Higgs Centre for Theoretical Physics, Edinburgh, United Kingdom

⁶Kavli Institute for Astronomy and Astrophysics, Beijing 100871, China

⁷William H. Miller III Department of Physics and Astronomy,
Johns Hopkins University, Baltimore, Maryland 21218, USA

⁸Astrophysics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, United Kingdom

⁹School of Mathematical Sciences and STAG Research Centre,
University of Southampton, Southampton, SO17 1BJ, United Kingdom

¹⁰SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste
IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

¹¹Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), D-14476 Potsdam, Germany

¹²Laboratoire Univers et Théories, CNRS, Observatoire de Paris,
Université PSL, Université Paris Cité, 92190 Meudon, France

¹³Theoretisch-Physikalisches Institut, Friedrich-
Schiller-Universität Jena, 07743, Jena, Germany

¹⁴Dipartimento di Fisica “G. Occhialini”, Università degli Studi
di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy

¹⁵INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy

¹⁶Institute for Mathematics, Astrophysics and Particle Physics,
Radboud University, 6525 AJ Nijmegen, The Netherlands

¹⁷The University of Arizona, Department of Astronomy, Tucson, AZ 85721, USA

¹⁸CENTRA, Departamento de Física, Instituto Superior Técnico - IST,
Universidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal

¹⁹Princeton Gravity Initiative, Princeton University, Princeton, New Jersey 08544, USA

²⁰Institut für Theoretische Physik, Universität Leipzig,
Brüderstraße 16, 04103 Leipzig, Germany.

²¹Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, CEP 22290-180, Brazil.

²²Department of Astronomy, Cornell University, Ithaca, New York 14853, USA

²³Department of Physics, National University of Singapore, Singapore 117551

²⁴Department of Mathematics, National University of Singapore, Singapore 119076

²⁵School of Mathematical Sciences, Queen Mary University
of London, Mile End Rd, London E1 4NS, United Kingdom

²⁶Departament de Física, Universitat de les Illes Balears, IAC3
- IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain
²⁷Université Libre de Bruxelles, International Solvay
Institutes, CP 231, B-1050 Brussels, Belgium

²⁸Université Libre de Bruxelles, Gravitational Wave Centre,
International Solvay Institutes, CP 231, B-1050 Brussels, Belgium

²⁹Center of Gravitational Physics, Weinberg Institute,
University of Texas at Austin, TX, 78712, USA
<sup>30</sup>GReCO, Institut d'Astrophysique de Paris, UMR 7095, CNRS,
Sorbonne Universit'e, 98bis boulevard Arago, 75014 Paris, France

<sup>31</sup>Department of Physics and Illinois Center for Advanced Studies of the
Universe, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA

<sup>32</sup>Department of Mathematics and Center for Scientific Computing &
Data Science Research, University of Massachusetts, Dartmouth, MA 02747

<sup>33</sup>Department of Physics and Astronomy, Vanderbilt University, Nashville TN 37235

<sup>34</sup>Instituto de Física Teórica IFT-UAM/CSIC,
Universidad Autónoma de Madrid, 28049 Madrid, Spain

<sup>35</sup>School of Physics and Astronomy & Institute for Gravitational Wave
Astronomy, University of Birmingham, Birmingham, B15 2TT, United Kingdom

<sup>36</sup>School of Mathematical Sciences, University of Nottingham
University Park, Nottingham NG7 2RD, United Kingdom

<sup>37</sup>Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands

<sup>38</sup>School of Physics and Astronomy, Cardiff University,
Queens Buildings, Cardiff, CF24 3AA, United Kingdom

<sup>39</sup>Departament de Física, Universitat de les Illes Balears, IAC3
– IEEC, Crta. Valldemossa km 7.5, E-07122 Palma, Spain

<sup>40</sup>University of Guelph, Guelph, Ontario N1G 2W1, Canada

<sup>41</sup>Institute for Theoretical Physics, Utrecht University,
Princetonplein 5, 3584CC Utrecht, The Netherlands

<sup>42</sup>Institute for Advanced Computational Science, Stony
Brook University, Stony Brook, NY 11794 USA

<sup>43</sup>Department of Physics and MIT Kavli Institute, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA

<sup>44</sup>Institut de Ciències de l'Espai (ICE, CSIC), Campus UAB,
Carrer de Can Magrans s/n, 08193 Cerdanyola del Vallès, Spain

<sup>45</sup>NASA Marshall Space Flight Center, Huntsville, Alabama 35811, USA

<sup>46</sup>Department of Physics and Center for Computational
Research, University of Rhode Island, Kingston, RI 02881

<sup>47</sup>Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca NY 14853

<sup>48</sup>Max-Planck-Institut für Astrophysik, Karl-
Schwarzschild-Straße 1, 85741 Garching, Germany

<sup>49</sup>Institute for Theoretical Physics, KU Leuven,
Celestijnenlaan 200D, B-3001 Leuven, Belgium

<sup>50</sup>Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

<sup>51</sup>Theoretical Particle Physics and Cosmology Group, Physics Department,
Kings College London, Strand, London WC2R 2LS, United Kingdom

<sup>52</sup>Center for Theoretical Astrophysics, Los Alamos
National Laboratory, Los Alamos, NM 87545, USA

<sup>53</sup>Center for Computational Relativity and Gravitation, School
of Mathematical Sciences, Rochester Institute of Technology,
85 Lomb Memorial Drive, Rochester, New York 14623, USA

<sup>54</sup>Nicholas and Lee Begovich Center for Gravitational-Wave Physics and Astronomy,
California State University, Fullerton, Fullerton, California, 92831, USA

<sup>55</sup>Astronomical Institute of the Czech Academy of Sciences,
Boční II 1401/1a, CZ-141 00 Prague, Czech Republic
⁵⁶Niels Bohr International Academy, Niels Bohr

Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark

⁵⁷School of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia

⁵⁸OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery, Australia

⁵⁹Gran Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy

⁶⁰INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy

⁶¹Department of Physics and Astronomy, Brigham Young University, Provo, UT 84602, USA

⁶²INFN Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

⁶³Institut des Hautes Etudes Scientifiques, 91440 Bures-sur-Yvette, France

⁶⁴Department of Physics, University of Virginia,
P.O. Box 400714, Charlottesville, VA 22904-4714, USA

⁶⁵Department of Physics, Faculty of Mechanical Engineering, Czech Technical
University in Prague, Technická 1902/4, Prague 6, 16607, Czech republic

⁶⁶Department of Physics, Pasadena City College, Pasadena, California 91106, USA

⁶⁷Center for Gravitational Physics and Quantum Information (CGPQI), Yukawa
Institute for Theoretical Physics (YITP), Kyoto University, 606-8502, Kyoto, Japan

⁶⁸The Hakubi Center for Advanced Research, Kyoto University,
Yoshida Ushinomiyacho, Sakyo-ku, Kyoto 606-8501, Japan

⁶⁹RIKEN iTHEMS, Wako, Saitama, 351-0198, Japan

⁷⁰Dipartimento di Fisica, Sapienza Università di Roma & INFN
Sezione di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy

⁷¹Department of Physics, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece

⁷²The University of Arizona, Departments of Astronomy and Physics, Tucson, AZ 85721, USA

⁷³Departamento de Astronomía y Astrofísica, Universitat de
València, Dr. Moliner 50, 46100, Burjassot (València), Spain

⁷⁴Institut d'Estudis Espacials de Catalunya (IEEC), Edifici Nexus,
Carrer del Gran Capità 2-4, despatx 201, 08034 Barcelona, Spain

⁷⁵University of Utah, Salt Lake City, UT 84112, USA

⁷⁶Laboratoire des 2 Infinis - Toulouse (L2IT-IN2P3), Université
de Toulouse, CNRS, UPS, F-31062 Toulouse Cedex 9, France

⁷⁷Theoretical Astrophysics Group, California Institute
of Technology, Pasadena, CA 91125, U.S.A.

⁷⁸Research Center for Astronomy and Applied
Mathematics, Academy of Athens, Athens 11527, Greece

⁷⁹Dipartimento di Fisica “E. Pancini”, Università di Napoli “Federico II” and INFN, Sezione
di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy.

⁸⁰School of Sciences & Humanities, Wentworth
Institute of Technology, Boston, MA 02155, USA

⁸¹Institute of Theoretical Physics, Faculty of Mathematics and
Physics, Charles University, CZ-180 00 Prague, Czech Republic

⁸²Centre for Research and Development in Mathematics and Applications,
Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal

⁸³Chennai Mathematical Institute, Siruseri, 603103, India

⁸⁴Nottingham Centre of Gravity, University of Nottingham,
University Park, Nottingham, NG7 2RD, United Kingdom

⁸⁵School of Mathematical Sciences, University of Nottingham,
University Park, Nottingham, NG7 2RD, United Kingdom
<sup>86</sup>Observatório Nacional, Rua General José Cristino 77,
São Cristóvão, 20921-400 Rio de Janeiro, RJ, Brazil

<sup>87</sup>SUPA, University of Glasgow, Glasgow G128QQ, United Kingdom

<sup>88</sup>Consortium for Fundamental Physics, School of Mathematics and Statistics, University of
Sheffield, Hicks Building, Hounsfield Road, Sheffield S3 7RH, United Kingdom

<sup>89</sup>INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807 Merate, Italy

<sup>90</sup>Department of Physics, Indian Institute of Science, Bangalore 560012, India

<sup>91</sup>Shanghai Astronomical Observatory, CAS, 200030, Shanghai, China

<sup>92</sup>University of Illinois

<sup>93</sup>Data Science and Learning Division, Argonne
National Laboratory, Lemont, Illinois 60439, USA

<sup>94</sup>Department of Computer Science, The University of Chicago, Chicago, Illinois 60637, USA

<sup>95</sup>Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, Illinois 61801, USA

<sup>96</sup>Department of Physics, University of Zürich,
Winterthurerstrasse 190, 8057 Zürich, Switzerland

<sup>97</sup>Center for Space Sciences and Technology, University of Maryland
Baltimore County, 1000 Hilltop Circle Baltimore, MD 21250, USA

<sup>98</sup>Gravitational Astrophysics Lab, NASA Goddard
Space Flight Center, Greenbelt, MD 20771, USA

<sup>99</sup>Center for Research and Exploration in Space Science and Technology,
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>100</sup>Computational Physics and Methods Group, Los
Alamos National Lab, Los Alamos, NM 87544, USA

<sup>101</sup>Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway

<sup>102</sup>Institute for Theoretical Physics, Goethe University, 60438 Frankfurt am Main, Germany

<sup>103</sup>Department of Physics and Astronomy, University
of Mississippi, University, MS 38677, USA

<sup>104</sup>Department of Applied Mathematics and Theoretical Physics, University
of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

<sup>105</sup>School of Mathematics and Statistics, University
of Canterbury, Christchurch 8041, New Zealand

<sup>106</sup>CEICO, Institute of Physics of the Czech Academy
of Sciences, Na Slovance 2, 182 21 Praha 8, Czechia

(Dated: 24 October 2024)
# Contents

<table><tr><td>Notation</td><td>9</td></tr><tr><td>List of acronyms</td><td>9</td></tr><tr><td>1 Introduction</td><td>11</td></tr><tr><td>2 LISA Sources</td><td>14</td></tr><tr><td>2.1 Massive black hole binaries (MBHBs)</td><td>14</td></tr><tr><td>2.2 Extreme mass-ratio inspirals (EMRIs)</td><td>18</td></tr><tr><td>2.3 Intermediate mass-ratio inspirals (IMRIs)</td><td>21</td></tr><tr><td>2.4 Galactic Binaries (GBs)</td><td>23</td></tr><tr><td>2.5 Stellar origin black hole binaries (SOBHBs)</td><td>27</td></tr><tr><td>2.6 Cosmic strings</td><td>29</td></tr><tr><td>2.7 Beyond GR and beyond Standard Model sources</td><td>33</td></tr><tr><td>3 Modelling requirements from data analysis</td><td>38</td></tr><tr><td>3.1 Data analysis for LISA</td><td>38</td></tr><tr><td>3.2 Accuracy requirements</td><td>41</td></tr><tr><td>3.3 Efficiency considerations</td><td>43</td></tr><tr><td>3.4 Interface and data format requirements</td><td>45</td></tr><tr><td>4 Modelling approaches for compact binaries</td><td>47</td></tr><tr><td>4.1 Numerical relativity</td><td>48</td></tr><tr><td>4.2 Weak field approximations (post-Newtonian/post-Minkowskian)</td><td>62</td></tr><tr><td>4.3 Small-mass-ratio approximation (gravitational self-force)</td><td>74</td></tr><tr><td>4.4 Perturbation theory for post-merger waveforms (quasi-normal modes)</td><td>85</td></tr><tr><td>4.5 Effective-one-body waveform models</td><td>90</td></tr><tr><td>4.6 Phenomenological waveform models</td><td>102</td></tr><tr><td>5 Waveform generation acceleration</td><td>111</td></tr><tr><td>5.1 Computational techniques</td><td>111</td></tr><tr><td>5.2 Hardware accelerators / configurations</td><td>115</td></tr><tr><td>6 Modelling for beyond GR, beyond Standard Model, and cosmic string sources</td><td>118</td></tr><tr><td>6.1 Numerical relativity</td><td>118</td></tr><tr><td>6.2 Weak-field approximations</td><td>122</td></tr><tr><td>6.3 Perturbation theory for post-merger waveforms</td><td>124</td></tr><tr><td>6.4 Small mass-ratio approximation</td><td>126</td></tr><tr><td>6.5 Effective-one-body waveform models</td><td>129</td></tr><tr><td>6.6 Phenomenological waveform models</td><td>130</td></tr><tr><td>6.7 Modeling cosmic strings</td><td>131</td></tr></table>
<table><tr><td>7 Afterword</td><td>135</td></tr><tr><td>A Descriptions of NR codes</td><td>137</td></tr><tr><td>Acknowledgements</td><td>140</td></tr><tr><td>References</td><td>143</td></tr></table>
## Notation

<table><caption>**Table 1.** Frequently used symbols throughout this whitepaper, *G* and *c* are set to 1 throughout.</caption><tr><td>G</td><td>gravitational constant</td></tr><tr><td>c</td><td>speed of light</td></tr><tr><td>m<sub>1</sub></td><td>mass of the primary</td><td>m<sub>1</sub> &gt; m<sub>2</sub></td></tr><tr><td>m<sub>2</sub></td><td>mass of secondary</td><td>m<sub>2</sub> &lt; m<sub>1</sub></td></tr><tr><td>M</td><td>total mass</td><td>m<sub>1</sub> + m<sub>2</sub></td></tr><tr><td>&mu;</td><td>reduced mass</td><td>m<sub>1</sub>m<sub>2</sub>/M</td></tr><tr><td>q</td><td>large mass ratio</td><td>m<sub>1</sub>/m<sub>2</sub></td></tr><tr><td>&epsilon;</td><td>small mass ratio</td><td>m<sub>2</sub>/m<sub>1</sub></td></tr><tr><td>&nu;</td><td>symmetric mass ratio</td><td>m<sub>1</sub>m<sub>2</sub>/M<sup>2</sup></td></tr><tr><td>S<sub>i</sub></td><td>spin angular momentum of body i</td><td></td></tr><tr><td>a<sub>i</sub></td><td>Kerr spin parameter of body i</td><td>S<sub>i</sub>/(m<sub>i</sub>c)</td></tr><tr><td>&chi;<sub>i</sub></td><td>dimensionless spin of body i</td><td>a<sub>i</sub>c<sup>2</sup>/(Gm<sub>i</sub>) &equiv; S<sub>i</sub>c/(Gm<sub>i</sub><sup>2</sup>)</td></tr><tr><td>e</td><td>orbital eccentricity</td><td></td></tr><tr><td>&Omega;</td><td>orbital frequency</td><td></td></tr><tr><td>&mu;<sub>S</sub></td><td>string tension</td><td></td></tr></table>

## List of acronyms

<table><tr><td>AAK</td><td>Augmented Analytic Kludge</td></tr><tr><td>AGN</td><td>active galactic nuclei</td></tr><tr><td>BBH</td><td>binary black hole</td></tr><tr><td>BH</td><td>black hole</td></tr><tr><td>BHNS</td><td>binary black hole-neutron star</td></tr><tr><td>BHPT</td><td>black hole perturbation theory</td></tr><tr><td>BNS</td><td>binary neutron star</td></tr><tr><td>BWD</td><td>binary white dwarf</td></tr><tr><td>CHE</td><td>close hyperbolic encounters</td></tr><tr><td>CO</td><td>compact object</td></tr><tr><td>CPU</td><td>central processing unit</td></tr><tr><td>CS</td><td>cosmic string</td></tr><tr><td>EFT</td><td>effective field theory</td></tr><tr><td>EMRI</td><td>extreme mass-ratio inspiral</td></tr><tr><td>EOB</td><td>effective one-body</td></tr><tr><td>FD</td><td>frequency domain</td></tr><tr><td>FEW</td><td>Fast EMRI Waveforms</td></tr><tr><td>GB</td><td>Galactic binary</td></tr><tr><td>GPU</td><td>graphical processing unit</td></tr><tr><td>GR</td><td>general relativity</td></tr><tr><td>GSF</td><td>gravitational self-force</td></tr><tr><td>GW</td><td>gravitational wave</td></tr><tr><td>IMBH</td><td>intermediate mass black hole</td></tr><tr><td>IMRI</td><td>intermediate mass-ratio inspiral</td></tr><tr><td>ISCO</td><td>innermost stable circular orbit</td></tr><tr><td>LISA</td><td>Laser Interferometer Space Antenna</td></tr><tr><td>LVK</td><td>LIGO-Virgo-KAGRA Collaboration</td></tr><tr><td>MBH</td><td>massive black hole</td></tr><tr><td>MBHB</td><td>massive black hole binary</td></tr><tr><td>NR</td><td>numerical relativity</td></tr><tr><td>NS</td><td>neutron star</td></tr></table>
**ODE** ordinary differential equation

**PBH** primordial black holes

**PDE** partial differential equation

**PM** post-Minkowskian

**PN** post-Newtonian

**QCD** quantum chromodynamics

**QNM** quasi-normal mode

**ROM** reduced-order model

**SDSS** Sloan Digital Sky Survey

**SNR** signal-to-noise ratio

**SOBH** stellar-origin black hole

**SOBHB** stellar-origin black hole binary

**SPA** stationary phase approximation

**TD** time domain

**UCXB** ultra-compact X-ray binary

**WDBH** white dwarf-black hole

**WDNS** white dwarf-neutron star

**XMRI** extremely-large mass-ratio inspiral
# 1. Introduction

Waveforms, theoretical predictions for gravitational wave (GW) signals, play a vital role in GW astronomy. Most GW signals are buried deep in instrumental noise. By using waveforms as matched filters, such signals can be detected. Once a GW is found the properties of its source can be inferred by further comparing the system to theoretical waveforms. The Laser Interferometer Space Antenna (LISA) will open a window on a new frequency band of GWs in the mHz regime. In this band we will encounter a slew of new sources of GWs ranging from local Galactic white-dwarf binaries to distant mergers of massive black holes (MBHs), and many sources in between. Realizing the science potential of LISA detection and measurement of GW signals will require new waveform models, which will need to cover a much wider range of sources while being significantly more accurate than the models used currently in ground based GW observations by LIGO, Virgo, and KAGRA [1, 2].

A large class of GW sources expected for LISA feature the inspiral of a binary of compact objects decaying under the emission of GWs, often leading to a merger at the end. Producing theoretical waveforms for such events requires solving the relativistic two body problem. This is a notoriously difficult problem with no known closed form solution for radiating binaries. There are three main approaches to obtain approximate solutions to the relativistic two-body problem from first principles. First, numerical relativity (NR) takes the non-linear partial differential equations (PDEs) of general relativity (GR), puts them on a grid, and evolves them on a supercomputer. Second, one can make a weak field approximation and obtaining a perturbative analytical solution leading to post-Newtonian (PN) or post-Minkowskian (PM) theory (depending on whether a simultaneous slow motion approximation is introduced). Finally, gravitational self-force (GSF) theory expands the two-body dynamics in powers of the (small) ratio of the masses of the two bodies. Studies of homogeneous perturbations of isolated black holes (BHs) and associated quasi-normal mode (QNM) calculations are used to model the post-merger behaviour.

Each of these four approaches has its own natural domain of applicability. Due to its computational cost NR is limited to a relatively low number of orbits before merger and to systems with comparable masses. The weak-field approximation of PN theory limits its application to the early inspiral, while GSF theory is at its best when dealing with binaries with relatively disparate mass-ratios. The post-merger phase of binaries containing a BH can be modelled with black hole perturbation theory (BHPT). As illustrated in Fig. 1, the natural domains of applicability of the four first-principle approaches to the relativistic two-body problem are largely complementary, and building waveforms for LISA will require input from all four — sometimes for modelling a single source. Effective waveform models like effective one-body (EOB) and the Phenom-family combine inputs from the different approaches to provide waveforms that can straddle these domains.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A sketch of the natural domains of applicability of the four main approaches to solving the relativistic two-body problem. The approaches are largely complementary and building waveforms for LISA will require input from all four. The solid lines shown are illustrative of the reach of each approach. The precise reach of each region depends on the source type and the accuracy requirements of the model. In addition, the Effective-One-Body and Phenomenological models absorb information from the four main approaches to produce more global models that can compute waveforms sourced by binaries across large portions of the parameter space.</figcaption></figure>

The LISA Waveform Working Group (WavWG) serves as an interface between the LISA Science Group — tasked with realizing LISA’s scientific mission — and the wider scientific community developing waveform models and studying the relativistic two-body problem. It serves to prioritize waveform development, inform the wider community of LISA’s waveform modelling needs, and as a recruiting pool for waveform related tasks and projects in the LISA Science Group.

In this white paper, the WavWG discusses the current status of waveform modelling and what further development is needed in realize LISA’s science. It is organized in five main sections. In §2 we take a brief inventory of the sources LISA is expected to observe along with their expected parameters. This sets the primary goals for waveform model accuracy and parameter space coverage. §3 discusses what requirements LISA data analysis puts on waveform models in terms of accuracy, efficiency, and format. The main approaches to modelling waveforms from compact binaries are described in §4, discussing both their status to date, and challenges to be overcome to meet LISA’s waveform requirements. §5 discusses methods for accelerating the production of waveform models.
Finally, §6 covers the modeling for beyond GR, beyond Standard Model, and cosmic strings sources. We note that while stochastic signals are of considerable interest, this paper does not address them as the focus here is on individual detectable sources. For the reader's convenience, we have provided a table of notation used throughout the paper in Table 1 and a list of commonly used acronyms on the page proceeding this introduction.

The writing of this white paper was coordinated by the co-chairs of the LISA Waveform Working Group: Maarten van de Meent, Deirdre Shoemaker, Niels Warburton, and Helvi Witek. Additional coordination was provided by the co-chairs of the LISA Waveform Work Package: Leor Barack, Anna Heffernan, and Harald Pfeiffer.

The coordinators and contributors to individual sections of the white paper are listed at the start of each section.
## 2. LISA Sources

LISA will be sensitive to a wide range of GW sources in the mHz regime. In this section, we take an inventory of anticipated sources that need to be modelled to extract useful scientific results and are well-enough understood at this stage to be modelled in some detail. A more extensive survey of LISA's expected sources and the science that can be done with them can be found in the whitepaper produced by the LISA Astrophysics Working Group [3]. The short summaries here mainly focus on the expected parameter ranges to provide a context for the modelling approaches discussed in the rest of the white paper. Conventional long-lived astronomical sources include massive black hole binaries (MBHBs), extreme mass-ratio inspirals (EMRIs), intermediate mass-ratio inspirals (IMRIs), Galactic binaries (GBs) and stellar-origin black hole binaries (SOBHBs). These sections (2.1–2.5) present LISA sources in GR and the Standard Model. Theoretically interesting but empirically speculative sources include cosmic strings (CSs) and sources based on physics beyond GR and beyond the Standard Model are discussed in Secs. 2.6 and 2.7. These sources, and resulting science, will also be applicable to other mHz gravitational wave detectors [4, 5].

A key metric for the detectability of sources is their signal-to-noise ratio (SNR). This is defined via the standard formula

$$ \text{SNR}^2 = 4 \int \frac{|h(f)|^2}{S_n(f)} df, \qquad (1) $$

where $h(f)$ is the GW strain and $S_n(f)$ the noise power spectral density of the detector [6, 7]. Unless otherwise stated this is computed for an observation period of up to 4 years.

### 2.1. Massive black hole binaries (MBHBs)

Coordinator: Enrico Barausse

Contributors: M. Bonetti, J. Garcia-Bellido

#### 2.1.1. Description
In the local Universe, MBHs are observed in the centers of virtually all massive galaxies [8, 9], as well as in many low mass dwarf galaxies [10–12]. Their cosmic evolution is inextricably intertwined with that of their host galaxies. The latter provide MBHs with gas out of which they grow by accretion (thereby shining as quasars or active galactic nuclei (AGNs)), and MBHs influence the evolution of their hosts by injecting energy into their surroundings (AGN feedback) [13–15]. This co-evolution of MBHs and galaxies is reflected in the scaling relations between MBH mass and galactic properties [16–21], although the detailed physical processes leading to their emergence are still not fully understood.
Present day galaxies/dark matter haloes are believed to have formed hierarchically from
the merger of smaller systems. Likewise, “seed” MBHs are expected to inhabit at least
a fraction of the high redshift progenitors of today’s galaxies (see [22] for a review).
These seed BHs are expected to grow by a combination of gas accretion and mergers
[23–29]. Their coalescences, as well as those of the later generations of MBHs that they
give rise to, are indeed a prime target for LISA [30–34].

By detecting MBH mergers, LISA could confirm or disprove this hierarchical scenario
for the evolution of MBHs. Indeed, considerable uncertainties exist about the timescale
on which MBHs form a bound binary and eventually coalesce, after two galaxies have
merged. The main uncertainty is whether MBH pairs can efficiently make their way
from separations of hundreds of pc down to the sub-pc scales on which gravitational
waves dominate the dynamics [33, 35–47]. First dynamical friction and then stellar
and/or gaseous hardening are the main drivers of the evolution of MBHs on their path
to coalescence. Triple or quadruple MBH interactions, which naturally arise in models
where the evolution of MBH pairs/binaries is slow or even stalls [48–50], should in any
case ensure LISA detection rates of at least a few per year [33, 48]. Conversely, the
observed detection rate will inform us about the efficiency with which MBH binaries
come together and merge in nature [33], even though that information may be partly
degenerate with the number of MBH seeds formed at high redshift. In that respect,
it is worth pointing out the possibility that there may be massive primordial black
holes (PBH) formed during the radiation era that may act as seeds for the MBH
population [51, 52]. In that case, MBH mergers could occur at much higher redshift.

2.1.2. Expected source parameters Further information on the underlying astrophysics can be gained by measuring the parameters of MBH binaries. Measurements of the component masses will allow for constraining the initial mass function of the seed black holes at high redshift [30–33, 53]. Indeed, the latter may form as relatively low-mass black holes (~100–1000$M_{\odot}$), e.g. from the collapse of the first generation of stars [54], or they may be born already as higher mass seeds of ~$10^4 – 10^5 M_{\odot}$ (e.g. from the collapse of massive quasi-stars [55], runaway instabilities in stellar clusters [56, 57], instabilities of protogalactic disks [23], etc.). Spin measurements, together with mass measurements, will clarify the properties of the accretion process and its importance in the evolution of MBHs relative to mergers [58–60]. Measurements of distance will allow for placing the detected system at the right epoch in the history of the universe. Further help in this respect will naturally be provided by sufficiently accurate/precise sky localization, which might allow for the identification of an electromagnetic counterpart; and, therefore, a direct measurement of photometric or spectroscopic redshift [61, 62]. In the presence of independent redshift and distance measurements, one may even attempt to construct a GW-based cosmography, potentially measuring the Hubble constant or the density of matter/dark energy in our universe [61]. Finally, measurements of the eccentricity and mass ratio will provide insight on the mechanism driving MBH pairs
to separations at which GW emission is enough to lead to a merger within a Hubble
time [33, 48].

In more detail, current expectations on the parameters of MBHBs and their seeds from
astrophysical models are generally as follows:

* Source frame total masses range from a few $10^2$ to a few $10^8 M_\odot$, with possible peaks around $10^3$ and $10^5 M_\odot$ for respectively light and heavy seed formation in some models. See, e.g., Fig. 2 of Ref. [33].

* Comparable mass binaries with mass ratio $\epsilon \equiv m_2/m_1$ between 0.1 and 1 are expected to make up the bulk of LISA detections, although tails extending down to $\epsilon \sim$ a few $\times 10^{-3}$ may be present (cf. Fig. 3 of [33]).

* Large uncertainties generally affect the prediction for spin magnitudes and orientations. Moderate to large spins $\gtrsim 0.7$ are expected in models calibrated to electromagnetic measurements (i.e. iron $K\alpha$) of MBH spins in local AGNs [60], but the spins of the first generations of MBH mergers are relatively unknown. Similarly, spins may be approximately parallel at the peak of the star formation $z \sim 2-3$, where gas torques may align them with the orbital angular momentum to within 10-30 degrees [63], but the situation may be different at higher and lower redshift.

* Eccentricity may be very significant and evolving in the LISA window ($e > 0.9$ upon entrance into the LISA band and up to 0.1 at plunge) if triple/quadruple MBH interactions are more efficient than gas and stellar interactions to drive the binary's evolution at $\sim$pc scales [48].

* Event rates and SNRs can vary significantly according to the underpinning astrophysical model [32, 33, 53]: in light seed scenarios rates may be as low as a few per year, especially if feedback from supernova explosions is included [33], which results in most events having small total masses and thus low SNRs (below 100). More generally, the number of light seed detections may be threshold limited [33, 48]. This is again especially true if supernova feedback is included, as most events have SNRs of a few tens in that case. In heavy seed scenarios, instead, irrespective of whether supernova feedback is included or not, rates and SNRs tend to be larger, i.e. roughly a few tens of detections in 4 yr and SNRs up to several thousands.

* Environmental effects from gas and stars generally have negligible impact in the LISA observation window [64], except under certain circumstances [65, 66]. In addition, they might leave a recognizable imprint in the eccentricity distribution of detected events [67]. The same applies to MBH triple/quadruple interactions, which are expected to cause the aforementioned very high “relic” eccentricity when the binaries enter the LISA band [48].

* When delays between galaxy and MBH mergers are accounted for, the highest number of detections is expected around $z \sim 2-4$ [33].
In addition to the large parameter space described above, waveform templates for MBHBs will also need to accurately represent the waveform for many hundreds or thousands of cycles before merger — see, e.g., Fig. 16 of Ref. [62].

<table><caption>**Table 3.** Summary of the anticipated parameters for LISA MBHB sources. The total mass of a massive binary in the source frame ranges from 10<sup>2</sup> – 10<sup>8</sup><i>M</i><sub>☉</sub> while moving to the detector frame picks up a factor of ~10 giving the range 10<sup>3</sup> – 10<sup>9</sup><i>M</i><sub>☉</sub>. This is squeezed on both ends by another factor of 10 (SNR and LISA band constrain the low and high limits respectively) giving a possible range of 10<sup>4</sup> – 10<sup>8</sup><i>M</i><sub>☉</sub> with the majority of sources expected between 10<sup>5</sup> – 10<sup>7</sup><i>M</i><sub>☉</sub> [33]. The mass ratio, <i>q</i> is expected to range from 1–1000, again with the majority of sources expected to be in the 1–10 range [33].</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Astrophysically relevant range</th></tr></thead><tbody><tr><td>Total mass in the detector frame</td><td><i>M</i></td><td>10<sup>5</sup> – 10<sup>7</sup><i>M</i><sub>☉</sub></td></tr><tr><td>Mass ratio (&gt; 1)</td><td><i>q</i></td><td>1 – 10</td></tr><tr><td>Dimensionless spin</td><td>max|<i>&chi;</i><sub>i</sub>|</td><td>0 – 0.998</td></tr><tr><td>Eccentricity entering LISA band</td><td><i>e</i><sub>init</sub></td><td>0 – 0.99</td></tr><tr><td>Eccentricity at last stable orbit</td><td><i>e</i><sub>merge</sub></td><td>&lt; 0.1</td></tr><tr><td>Signal to noise ratio</td><td>SNR</td><td>10 – 10<sup>4</sup></td></tr></tbody></table>
## 2.2. Extreme mass-ratio inspirals (EMRIs)

Coordinator: Huan Yang

Contributors: S. Akcay, P. Amaro Seoane, S. Bernuzzi, J. Gair, A. Heffernan, T. Hinderer, S. A. Hughes, S. Isoyama, G. Lukes-Gerakopoulos, A. Nagar, Z. Pan, Z. Sam, C. F. Sopuerta, V. Witzany

### 2.2.1. Description
EMRIs are binaries in which one member is substantially more massive than the other. The canonical EMRI is expected to be a stellar mass compact body (typically a black hole of $10-100 M_{\odot}$) captured by a BH of $10^5-10^7 M_{\odot}$ by multibody scattering processes in the core of a galaxy [68]. EMRIs can also form via the capture of other compact objects such as neutron stars or white dwarfs, though the small mass of these objects reduces their detectability with LISA [69]. A consequence of the large mass-ratio of EMRIs is the quite slow inspiral of the smaller body, executing $10^4-10^5$ of orbits as it moves through the black hole's strong field [70, 71]. Binaries involving a MBH orbited by a substellar object, such as a brown dwarf, are called extremely-large mass-ratio inspirals (XMRIs) and these are potential LISA sources if they form around the MBH in the centre of the Milky Way [72]. Measuring the gravitational waves generated by the strong-field orbits of E/XMRIs will make it possible to map the properties of massive black holes with great precision [73, 74].

### 2.2.2. Expected source parameters
The expected LISA detection rate of EMRIs depends on several astrophysical ingredients [73], including the mass/spin function of MBHs at different redshifts, the fraction of MBHs living in stellar cusps, stellar-mass compact objects (COs) capture rate per MBH, the characteristic mass of stellar-mass COs and their orbit eccentricities. The mass function of MBHs is usually modelled as a power law $dn/d \log M = n_0 (M/3 \times 10^6 M_{\odot})^{\pm0.3}$ [59, 75], with $n$ being the MBH number density and $n_0 \in (0.002, 0.005)$ Mpc$^{-3}$. MBHs can be extremely spinning with $a \approx 0.998$, slowly spinning with $\chi_1 \approx 0$ or of an extending spin distribution in the range $0 \le \chi_1 \le 1$, depending on their growth channels (accretion and/or mergers) [60, 76, 77]. The fraction of MBHs living in stellar cusps is determined by how frequently MBHs merge, when the stellar cusps are destroyed, and their regrowth time after mergers. COs in stellar cusps captured by MBHs generically reside in low angular momentum orbits with large eccentricities, which subsequently decay to the range of $0 < e < 0.2$ at the final plunge [73]. The capture rate depends on the cusp relaxation time and the density profile of the cusp [78]. Taking into account of all these uncertainties with semi-analytic models, Babak et al. [73] forecasted that there will be several to thousands of EMRIs detected by LISA per year, assuming the contributing COs are BHs with mass $10$ or $30 M_{\odot}$, and a detection threshold of SNR = 20. On the other hand, recent rate studies predicted that accretion disk-assisted EMRI formation may be more common for LISA detection [79–81], thanks to the high efficiency of disks in transporting stellar-mass
black holes towards galaxy centers. These disk-assisted EMRIs tend to have negligible eccentricity comparing to the EMRIs formed by gravitational capture.

For XMRIs their extremely slow evolution means they will spend millions of years in the LISA band. Although the event rate of an XMRI forming is quite low, their long duration in band means in the Galactic centre there could be a few dozen sources, some of which may be highly eccentric. Due to the short distance between the centre of the Galaxy and the solar system, the SNRs can range from $10$ to $10^4$ [72].

### 2.2.3. Science with EMRIs
Constraining model parameters of EMRIs from their GW signals benefits from the large number of cycles completed in the LISA band [73, 77, 82–84]: all intrinsic parameters of an EMRI (the redshifted masses, the MBH spin magnitude, and the CO orbit eccentricity at plunge) can be measured with fractional uncertainty $\sim 10^{-6} - 10^{-4}$; external parameters are mainly determined by the GW amplitude and its modulation with time, where the luminosity distance $D_L$, the source sky location $\Omega_s$ and the MBH spin direction $\Omega_a$ can be measured with fractional uncertainty $\sigma(\ln D_L) \approx \text{SNR}^{-1}$, median solid angle uncertainties $\sigma(\Omega_s) \approx 0.05(\text{SNR}/100)^{-5/2}\text{deg}^2$ and $\sigma(\Omega_a) \approx 30\sigma(\Omega_s)$, respectively. In light of the source location constraints in both radial and transverse directions, a fraction of low-redshift EMRIs are expected to be traced back to their host galaxies without the aid of any electromagnetic counterparts [77].

To identify the nature of the secondary in an EMRI, in the case of a light black hole or a neutron star, the only potentially detectable parameter beyond its mass will be a component of its spin [85, 86]. For less compact objects, such as white dwarfs or brown dwarfs, there is a possibility of gaining some information about their tidal quadrupoles [87]. However, all of the parameters of the secondary beyond its mass may be poorly constrained and partially or wholly degenerate with parameters controlling other sub-leading effects.

Detection of EMRIs will provide critical information with which to understand the growth history of MBHs and their environment. Different growth mechanisms, i.e., via mergers and/or via gas accretion, predict significantly different mass and spin distribution for MBHs [58–60, 88]. With a population of EMRI events, the co-evolution of MBH masses and spins over cosmic time can be measured to further reveal how they are formed and how they grew [58–60, 77]. N-body calculations have shown that EMRI rates sensitively depend on stellar distribution at galactic centers, so that the rate inferred from observation will also be able to constrain various distribution models [73]. Environmental effects, including the interaction with possible accretion disk around the MBH and close stellar objects near the EMRI system, may induce sizeable phase shifts to the EMRI waveform [64, 89–91]. The accretion channel has illustrated that MBHs with AGNs will account for a significant number of EMRIs [80]. Estimates of the resulting GW phase correction vary over many orders of magnitude ($0.1 - 10^4$ rad /
year), depending on the assumed disk model [92, 93]. Tidal resonances induced by the tidal field of nearby stars and stellar-mass black holes can generate a phase correction of $O(1)$ in a significant fraction of the phase space [94]. The main uncertainty comes from the population of nearby stellar-mass objects, which may come from the (disk)-migrating stellar-mass objects in the AGN or leftovers from the previous AGN life cycle (Fig. 6 of [95]), or simply from the mass segregation effect ([96, 97]).

By precisely mapping the background spacetime geometry of a MBH [98], EMRIs can be used to check whether the Kerr metric description is accurate and GR holds in the strong-field regime [99]. At the innermost stable circular orbit (ISCO) of a typical MBH, the curvature scale is higher than Solar-system curvatures but lower than curvatures of some well-observed compact-binaries [100]. Consequently, modified gravities satisfying existing observational constraints generically predict very small deviations from GR in the standard EMRI scenario. However, since EMRIs are expected to be highly sensitive to small departures from the standard Kerr black hole background paradigm [101, 102], EMRIs will allow us to further tightly constrain any deviations caused by an additional scalar or vector channel for radiation in the long-term dissipation of energy and angular momentum [103]. In fact, even in the case where the field of the primary is just a generic stationary axisymmetric GR vacuum field, we should be able to determine the matter multipoles of the primary [104, 105] and check possible violations of the no-hair theorem [106], e.g., a deviation from the Kerr solution in the quadrupole moment [107].

In addition, as is true for all extragalactic LISA sources, various cosmological dispersion and propagation effects can be tightly constrained from EMRI signals, especially if the EMRI host galaxy is identified [103].

<table><caption>**Table 4.** Summary of the anticipated parameters for LISA EMRI sources. Total masses of $10^5 – 10^7 M_{\odot}$ will be seen by LISA [75] with those greater than $10^7 M_{\odot}$ starting outside the LISA band [68]. Regardless of formation channel, the smaller compact body is expected to be a stellar mass black hole giving an expected mass ratio of $10^4 – 10^6$ [73, 79], or even $10^4 – 10^7$ if neutron star EMRIs are included. Using MBH of $3 \times 10^6 M_{\odot}$ as a basis, the “loss-cone” EMRIs are expected to have high eccentricities, capped at ~ 0.8 entering the LISA band [108], while the accretion channel tends to have smaller eccentricities [80]. Evolving a large sample of compact bodies from capture results in a flat eccentricity distribution for the last stable orbit in the range $0 < e < 0.2$ [73]. A small tail of outlying higher eccentricities is possible. There is also the possibility of GW bursts from unbound / hyperbolic systems [109, 110]. Ref. [73] predicts the loudest EMRIs to have SNR ~ 100 per year.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Astrophysically relevant range</th></tr></thead><tbody><tr><td>Total mass in the detector frame</td><td><i>M</i></td><td>10<sup>5</sup> – 10<sup>7</sup><i>M</i><sub>☉</sub></td></tr><tr><td>Mass ratio (&gt; 1)</td><td><i>q</i></td><td>10<sup>4</sup> – 10<sup>6</sup></td></tr><tr><td>Dimensionless spin</td><td>max|<i>&chi;</i><sub>i</sub>|</td><td>0 – 0.998</td></tr><tr><td>Eccentricity entering LISA band</td><td><i>e</i><sub>init</sub></td><td>0 – 0.8</td></tr><tr><td>Eccentricity at last stable orbit</td><td><i>e</i><sub>merge</sub></td><td>0 – 0.2</td></tr><tr><td>Signal to noise ratio</td><td>SNR</td><td>20 – 100</td></tr></tbody></table>
## 2.3. Intermediate mass-ratio inspirals (IMRIs)

Coordinator: Carlos F. Sopuerta

Contributors: M. Dhesi, M. Hannam, T. Hinderer, D. Neilsen, H. Pfeiffer, H. Sundar, N. Warburton

### 2.3.1. Description
An IMRI is a binary system with mass ratio in the range $10 \lesssim q \lesssim 10^4$ placing it between comparable mass binaries and EMRIs. IMRIs can come in two flavors, both of which contain an intermediate mass black hole (IMBH) with mass in the range $10^2-10^4 M_{\odot}$. There is little observational evidence for BHs in the IMBH mass range, mainly because their main formation channel is inside globular clusters, the interiors of which are very difficult to observe [111]. However, recent GW observations with ground-based detectors have found a population of BHs formed from binary black hole (BBH) mergers with a total mass greater than $10^2 M_{\odot}$ [112, 113]. The largest observed remnant (of GW190521) has a mass of $\sim 150 M_{\odot}$ [114, 115]. At the other end of the mass range, we have evidence of BHs with masses in the range $10^4-10^5 M_{\odot}$ from observations of dwarf galaxies out to redshift $\sim 2.4$ [116, 117]. These may correspond to an interpolation of the mass function observed for low MBH mass [11].

The two flavors in which IMRIs are expected to appear are: (i) Light IMRIs. A stellar-origin black hole (SOBH), or another sufficiently massive and compact object, inspiraling into an IMBH. For instance, dwarf galaxies or globular clusters may contain an IMBH which could capture an SOBH [118]. In this case, it is quite likely that the merger occurs outside the LISA frequency band, at higher frequencies. It is thus possible that an IMRI can be observed during its inspiral with LISA and have its merger seen by ground-based GW observatories (only for light IMBHs, with mass $\lesssim 10^3 M_{\odot}$ [119]). (ii) Heavy IMRIs. Dynamical friction can produce the orbital decay of globular clusters into a galactic nucleus allowing the formation of an IMBH-MBH binary system, with the MBH in the LISA mass range ($10^5-10^7 M_{\odot}$) [120, 121]. Another possible channel to form heavy IMRIs is the merger of a dwarf galaxy satellite with its main galaxy [122]. A third possibility is the formation of an IMBH in a galactic nuclei via hierarchical mergers [123] via migration traps in AGN [124]. Interestingly, GW190521 may have occurred in such an environment [125].

For more information on possible formation channels, and the uncertain event rates for these binaries, see the LISA Astrophysics Working Group White Paper [3] and [126-128], or the IMRI reviews [68, 126, 129]. The nature of the different formation channels for IMRIs (both light and heavy) may lead to distinct IMRI dynamics so that LISA detections may shed some light on what are the most viable formation scenarios.

### 2.3.2. Science with IMRIs
As we only have evidence for IMBHs at the two ends of the mass range that defines these objects, the search for IMRIs with LISA is especially
relevant with huge discovery potential. Event rates and parameter estimates for IMRIs can provide valuable information about the formation and growth of IMBHs in globular clusters, as well as details of stellar dynamics in those systems. This information is often difficult to glean from electromagnetic observations [111].

For heavy IMRIs the masses and spin of the MBH can be measured to within an relative error of a fraction of a percent [130]. Unlike EMRIs, where the spin on the secondary can be difficult to constrain [131, 132], the spin of the IMBH can be constrained to within ~10% [130]. IMRIs spend a long time orbiting in the strong field and due to their high SNR (relative to EMRIs) this makes then uniquely precise probes of GR and the Kerr hypothesis [133]. They can also be sensitive to modifications to Einstein's equations of GR, such as modified gravity, extra fields, and extra dimensions [134, 135].

Finally, *multiband* GW astronomy observations of light IMRIs consisting of observing the inspiral in the LISA detector and the merger in ground-based detectors [136–139] offers a range of increased science potential [128, 140, 141]. In particular, observing the inspiral with LISA can place tight constraints on the sky location and the merger time which can allow for targeted alerts for electromagnetic counterparts to the merger. Conversely, knowledge of the source parameters from ground-based observations of the merger can allow researchers to perform targeted searches for IMRIs in the archived LISA data [142]. Multiband observations can also enhance the potential for performing tests of GR by breaking degeneracies between some parameters in the waveform models [143].

<table><caption>**Table 5.** Summary of the anticipated parameters for LISA IMRI sources. By definition, light IMRIs' primary BH is from $10^2 M_{\odot}$ to $10^4 M_{\odot}$, i.e., an IMBH, while heavy IMRIs' primary BH is a MBH [3]. Again the definition of an IMRI has a mass ratio from comparable ($10^2$) up to $10^5$, above which the binary is defined as an EMRI. Eccentricity can be quite high for both light and heavy IMRIs entering the LISA band [126]. SNRs are taken from Table 9 of the Astrophysics Working Group White Paper [3]. Although this table represents the best knowledge for each individual parameter, further restrictions of the IMRI parameter space may be possible by combining SNR information with astrophysical knowledge.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Heavy IMRIs</th><th>Light IMRIs</th></tr></thead><tbody><tr><td>Binary mass</td><td><i>M</i></td><td>10<sup>4</sup> – 10<sup>7</sup><i>M</i><sub>☉</sub></td><td>10<sup>2</sup> – 10<sup>4</sup><i>M</i><sub>☉</sub></td></tr><tr><td>Mass ratio (&gt; 1)</td><td><i>q</i></td><td>10 – 10<sup>4</sup></td><td>10 – 10<sup>4</sup></td></tr><tr><td>Dimensionless spin</td><td>max|<i>χ</i><sub>i</sub>|</td><td>0 – 0.998</td><td>0 – 0.998</td></tr><tr><td>Eccentricity entering LISA band</td><td><i>e</i><sub>init</sub></td><td>0 – 0.9995</td><td>0 – 0.9995</td></tr><tr><td>Eccentricity at last stable orbit or leaving LISA band</td><td><i>e</i><sub>merge</sub></td><td>0 – 0.9</td><td>0 – 0.9</td></tr><tr><td>Signal to noise ratio</td><td>SNR</td><td>10 – 10<sup>2</sup></td><td>10 – 10<sup>3</sup></td></tr></tbody></table>
## 2.4. Galactic Binaries (GBs)

Coordinator: Milton Ruiz

Contributors: V. Korol, H. Lim, V. Paschalidis, S. L. Shapiro, A. Tsokaros

**2.4.1. Description** Our Galaxy is home to a variety of stellar binaries formed by white dwarfs, neutron stars, and black holes. Approximately a million years before their anticipated merger, these compact binaries transition into the millihertz GW frequency band accessible with LISA. Given such long evolution timescales, these systems will manifest as nearly-monochromatic sources of gravitational waves. Furthermore, these compact binaries will be detectable in large numbers by LISA, potentially emerging as the most numerous GW source within the millihertz band [144].

Galactic compact binaries represent a key source class for the LISA mission. Firstly, their existence as GW sources is confirmed, as several have already been identified and characterised with electromagnetic telescopes (see [145] for overview). Gravitational radiation from some of these known sources has already been measured indirectly by monitoring binary orbital contraction over extended periods [146–148]. Secondly, as we can anticipate their GW signatures before the mission, these binaries have been suggested as “verification sources” for LISA (e.g. [149]). Most notably, Galactic compact object binaries accessible to LISA promise a wealth of insights into stellar and binary evolution. This encompasses understanding the nature of compact objects, unravelling the physical processes that govern binary interactions, and exploring their role in the formation and evolution of the Milky Way – see the LISA Astrophysics white paper [3] for a review. Finally, it is noteworthy that a substantial portion of these binaries will not be individually resolvable by LISA, contributing to an unresolved stochastic foreground that will act as an additional source of noise for the instrument. Therefore, accurate characterization of this foreground is crucial to ensuring the precise characterization of other LISA sources.

Here we consider Galactic double compact objects, which we refer to as GBs, specifically those where at least one component is a compact object other than a black hole (binaries with two stellar origin black holes are treated in the next section). These can be categorised into various subtypes: binary white dwarfs (BWDs), with the subset of accreting BWDs, known as AM Canum Venaticorum or AM CVns in literature; binary neutron stars (BNSs); and mixed systems. The mixed systems encompass white dwarf-neutron stars (WDNSs), which can emerge as ultra-compact X-ray binaries (UCXBs) in electromagnetic radiation, white dwarf-black holes (WDBHs), and binary black hole-neutron stars (BHNSs).

**2.4.2. Masses, mass-ratios, eccentricities, and known LISA verification sources** Observing GB using electromagnetic observatories poses a significant challenge due to their
inherently small size and dimness (in some cases the entire absence of electromagnetic emissions). These features, when combined with selection effects and incompleteness of dedicated electromagnetic surveys, has limited our ability to know the true distributions of their parameters such as orbital separations, component masses, mass ratios, and eccentricities. Thus, much of our understanding of these binaries primarily hinges upon population synthesis studies. However, electromagnetic observations using Sloan Digital Sky Survey (SDSS) and Gaia have recently improved our understanding [144, 150, 151]

Binary population synthesis studies (e.g., Table 6 in [3]) show that for frequencies less than approximately 2 mHz, the BWDs form an unresolved foreground for LISA (see also [152–155]). Binaries with frequencies greater than approximately 2 mHz and/or closer than a few kpc do not overlap with this background [3, 152, 156] and are called “resolvable”. However, not all resolvable GBs are detectable, and unresolved GBs can also be detectable. Only binary systems with a significant SNR may be detected by LISA within a 4-year mission. Theoretical estimates suggest that among the hundreds of millions of binaries in the Galaxy, the number of resolvable binaries are about 6,000-10,000 BWDs, 100-300 WDNSs, 2-100 BNSs, 0-3 WDBHs and 0-20 BHNSs [3]. However, the number of detectable binaries above the noise are about 6000 BWDs, 100 WDNSs, 30 BNSs, 3 WDBHs and 3 BHNSs [157].

The chirp mass, a key distinctive parameter among compact binary classes, varies significantly across different binary systems. For instance, BWDs typically peak at a chirp mass of around $0.25M_{\odot}$ (e.g. [158]), while BNS tend to have a chirp mass around $1.2M_{\odot}$ [159]. Systems comprising BHs and neutron stars (NSs) exhibit even higher chirp masses. However, it is important to note that LISA’s limited ability to measure frequency derivatives for $f < 2$ mHz, and thus the chirp mass may introduce potential classification ambiguities. For example, a nearby BWD system may be misidentified as a more distant BNS. Another discriminative feature could be the eccentricity measurements [159, 160]. BWDs, typically formed via isolated binary evolution, are expected to be circularised due to recurrent mass transfer episodes. Contrastingly, BNS systems detectable in the LISA band might present measurable eccentricities due to natal kicks imparted by the supernova explosions that birth the NSs. Nevertheless, some rare eccentric BWD can form in globular clusters or via triple interactions (e.g. [161]).

There are approximately 30 GBs, identified through electromagnetic detection, that have been confirmed as resolvable verification targets for the LISA mission. Notably, 16 of these are expected to be detectable within merely three months of scientific operations, proving invaluable in the early stages of the mission by assisting in the validation of LISA’s performance relative to pre-launch expectations. These GBs are characterised by orbital periods ranging from ~300 to ~6000 seconds, equivalent to GW frequencies approximately between 0.5 and 6 mHz. Their total mass typically lies between 0.5 to $1.0M_{\odot}$, with mass ratios roughly between 0.1 and 0.7. Their estimated SNRs can exceed 100 after four years of integration [162].
Certain types of binaries can be sources of low and high frequencies, and hence can straddle both the LISA and ground-based detector bands. This is true for mixed binaries when unstable mass transfer leads to merger. Gravitational waves from the inspiral and merger of WDNSs may go from LISA to ground-based from potential oscillations of the NS after merger and/or its eventual collapse to a BH [163-165]. The same straddling of LISA/ground-based frequency bands holds for stellar-mass black hole binaries [141].

### 2.4.3. Modeling requirements & methods to improve/remove the Galactic background

GBs can be classified as detached or interacting (mass-transferring), depending on whether significant mass transfer occurs between the components. This classification is crucial for modelling their GW signals accurately. Although these systems are practically monochromatic to first order, the frequency does evolve with time. In particular, there is a difference in how the frequency evolves for detached binaries and interacting binaries. The orbital evolution of detached binaries is driven by the emission of gravitational waves, causing a gradual inward spiral of the binary (manifested as a negative frequency derivative in the GW signal). On the other hand, in accreting binaries, the redistribution of mass between the components results in an increase in their separation (leading to a positive frequency derivative in the GW signal) [163, 166]. This leads to different requirements for waveform modelling. Note also that the shortest period binaries can be tidally-locked and corotating [167]. For binaries near the Roche separation, tidal effects may have a non-negligible contribution to the phase evolution of the wave [168, 169]. Understanding the frequency evolution will help determine if mass transfer is taking place [170] and/or tidal interactions are at play. A large number of mass-transferring systems may constrain the physics of mass transfer and the efficiency of angular momentum removal from the disk/companion system and its reinjection back into the orbit [163].

Some of the GW sources that LISA will observe may be part of triple or higher-order multiple systems (e.g. [171, 172]). This includes sources in the Galactic disc that form through isolated triple evolution, as well as those in dense environments. When a BWD system — as the most common Galactic LISA source — is part of a larger system with an additional stellar or even a substellar object, the gravitational interaction with the substellar object introduces a modulation in the observed GW frequency due to the Doppler and related effects [173]. This is a consequence of the motion of the BWD around the centre of mass of the three-body system. As the orbit of the tertiary must be larger than that of the inner binary for the system to remain stable, the modulation timescale will be longer than the GW frequency produced by the inner BWD. However, to detect this modulation with LISA, it should also be shorter than the mission's lifetime [174-176].

Signals from many GBs will be below the BWD confusion background, which is generated both by detached and interaction/mass-transferring BWDs. Efforts have been made to remove this background in order to resolve more binaries (see e.g. [177-179]).
[180]). Notice that the modeling of this confusion background may provide insight about the BWD population in the Milky Way (and nearby satellite galaxies).

<table><caption>**Table 6.** Summary of the anticipated parameters for LISA GB sources. BWDs can have a chirp mass between 0.1 – 1<i>M</i><sub>☉</sub> [181], mass ratios of 1 – 10 (which can reach 100 for AM CVns systems) and spins ranging from seconds to hours [182, 183]. WDNSs and BNSs will have chirp masses in the range 0.4 – 1.2<i>M</i><sub>☉</sub> [184] and 1.1 – 1.6<i>M</i><sub>☉</sub> [185, 186] respectively. These binaries involving a pulsar have spins from seconds down to milliseconds [187] and near any eccentricity [159, 185]. As these sources are considered nearly monochromatic in the LISA band, their eccentricity is expected to evolve only a negligible amount inband.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>BWD</th><th>WDNS</th><th>BNS</th></tr></thead><tbody><tr><td>Total Chirp Mass</td><td><i>M</i></td><td>0.1 – 1<i>M</i><sub>☉</sub></td><td>0.4 – 1.2<i>M</i><sub>☉</sub></td><td>1.1 – 1.6<i>M</i><sub>☉</sub></td></tr><tr><td>Mass Ratio (&gt; 1)</td><td><i>q</i></td><td>1 – 10</td><td>1 – 5</td><td>1 – 1.6</td></tr><tr><td>Eccentricity in LISA band</td><td><i>e</i><sub>init</sub></td><td>0</td><td>0 – 1</td><td>0 – 1</td></tr><tr><td>Signal to Noise Ratio</td><td>SNR</td><td>&lt; 1000</td><td>&lt; 1000</td><td>&lt; 1000</td></tr></tbody></table>
## 2.5. Stellar origin black hole binaries (SOBHBs)

Coordinators: Antoine Klein and Ilya Mandel

Contributor: D. Gerosa

**2.5.1. Description** Stellar origin BH binaries are BBHs with component masses ranging from a few solar masses up to about ~50-100$M_{\odot}$, where models of pair-instability supernovae predict a mass gap to appear [188–190]. Those systems are in a mass range such that they can be observed by both LISA and by ground-based detectors, as they sweep through a few decades in frequency during the last stage of their inspiral. This broad coverage of the GW frequency spectrum makes it possible to probe the evolutionary history of such binaries [141, 191].

A circular BBH with two 45 $M_{\odot}$ components will merge in about 10 years from a GW frequency of ~10 mHz [141, 192]. It is thus possible to track its evolution from the LISA band to the ground-based detector band, taking advantage of LISA’s ability to better measure some of the system properties at low frequency and the improved measurements of parameters such as the remnant spin (estimated from the post-merger ringdown) from high-frequency data. Of particular interest for LISA is the sky localisation capability, which is generally quite limited for BBHs in current ground-based detector data [193]. The long duration of observations in the LISA band means that the detector will complete multiple orbits around the Sun, and thus effectively act as an instrument with a baseline of order the size of the orbit. The accuracy of sky localisation can be estimated as [191]

$$ \sigma_{\theta} \sim 0.025 \frac{0.01\text{Hz}}{f} \frac{8 \text{ AU}}{\text{SNR baseline}}, \quad (2) $$

where the GW frequency $f$ should be replaced by the detector bandwidth for a source that evolves out of the band during the observation, and SNR is the signal-to-noise ratio. For a 2 AU baseline, this would yield sub-degree localization, though the exact localization accuracy would depend on the source location in the orbital plane. Another example could involve the measurements of BH spins and spin-orbit misalignments, which store information about formation scenarios [123, 194–197]. While we may expect that high-frequency ground-based observations are best suited to measure mass ratios and spin-orbit couplings, which enter the waveform at higher orders in the orbital frequency, LIGO-Virgo-KAGRA Collaboration (LVK) observations to date have demonstrated the challenge of making such measurements precisely [113]. LISA could observe $\gtrsim 10^5$ orbital cycles and therefore possibly place more precise constraints.

Beyond individual sources, LISA could track changes in the source population as binaries evolve in frequency. For example, LVK observations indicate that the majority of black holes in merging BBHs have masses $\lesssim 40M_{\odot}$ [198], low to moderate spins [199–201], and most have circularized by the time of merger [202]. On the other hand, LISA could observe BBHs while they are still eccentric, which would likely indicate dynamical
formation [203, 204]. The appearance of sources at higher frequencies would on its own indicate high birth eccentricity [205].

Another possibility worth mentioning is the detectability of a third body, particularly a massive black hole (MBH), through its impact on the binary’s GW signature. Nuclear clusters around an MBH have been proposed as a possible BBH merger site, with possible assistance from gas in an active galactic nucleus [124, 206, 207]. Ref. [191] argued that even when the orbital period of a binary of mass $M_{\text{bin}}$ around an MBH of mass $M_{\text{MBH}} \gg M_{\text{bin}}$ is much longer than the observation duration $T_{\text{obs}}$, the orbital acceleration of the binary due to the MBH could still be detectable provided the GW frequency $f_{\text{GW}}$ exceeds

$$f_{\text{GW}} \gtrsim 0.02 \, \text{Hz} \left( \frac{M_{\text{MBH}}}{10^6 M_{\odot}} \right)^{-1} \left( \frac{a}{\text{pc}} \right)^2 \left( \frac{T_{\text{obs}}}{5\text{yr}} \right)^{-2} \left( \frac{\text{SNR}}{8} \right)^{-1}, \quad (3)$$

where $a$ is the distance from the binary to the MBH. A BBH merger in a massive globular cluster of a similar mass could carry a comparable signature.

### 2.5.2. Expected source parameters

Judging by evidence from ground-based observations, the majority of BBHs observed so far have chirp masses between 5 and $\sim 40M_{\odot}$ [198], though BH masses could extend down to the maximum neutron star mass and up to the IMBH mass range, especially if hierarchical mergers in dense dynamical environments fill the mass gap from pair-instability supernovae [123, 208, 209]. LISA will be particularly sensitive to more massive BBHs.

Most observed BBHs have mass ratios $q \lesssim 3.3$, which is consistent with the bulk of model predictions, but systems with more extreme mass ratios are possible, such as GW190814, which involved a $23 M_{\odot}$ BH and a $\sim 2.6 M_{\odot}$ companion [210]. As mentioned above, most observed BBHs have low to moderate companion spins, though it remains unclear whether this is a generic feature of stellar-mass BHs [211, 212]. Lastly, high eccentricities and generic spin-orbit misalignments could be a telltale sign of dynamical formation in dense stellar environments or in hierarchical triples (see Refs. [213, 214] for reviews of formation channels).

BBHs emit the bulk of their orbital energy above the LISA frequency band; therefore, moderate SNRs are expected except for fortuitously nearby sources [141]. The minimal SNR for detection of BBHs will be influenced by the technical challenges associated with searching for signals with a complicated morphology and many in-band cycles [215]. Furthermore, beyond an SNR threshold for detection, signals must also show evidence of frequency evolution in order for masses to be measurable, so that BBHs can be identified among the much larger population of signals from double white dwarfs [159].
<table><caption>**Table 7.** Summary of the anticipated parameters for LISA SOBHB sources. The most recent data release from LVK indicates a mass distribution of BBHs centred around 5–40*M<sub>☉</sub> [198], however their masses could theoretically extend down to neutron star mass (2*M<sub>☉</sub>) and up towards IMBH (100*M<sub>☉</sub>) [123, 190, 208, 209]. LVK observations also show a mass ratio and effective spin distribution between 1–1.33 and 0–0.3 respectively [198], with more extreme mass ratio binaries observed [210] and the possibility of higher spins theoretically possible [211, 212]. Sources that appear in the LISA band will fully circularise by the time they reach the LVK band, so LVK observations cannot constrain LISA band eccentricities, which depend on the still uncertain formation channels. If we consider a SNR threshold for detection of 12, then 98.7% of sources in an isotropic homogeneous Universe would have SNR < 50.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>range for majority of sources</th></tr></thead><tbody><tr><td>Chirp Mass</td><td><i>M</i></td><td>5–40</td></tr><tr><td>Mass Ratio (&gt; 1)</td><td><i>q</i></td><td>1–3.3</td></tr><tr><td>Dimensionless Spin</td><td>max|<i>&chi;</i><sub>i</sub>|</td><td>0–0.3</td></tr><tr><td>Eccentricity entering LISA band</td><td><i>e</i><sub>init</sub></td><td>0–1</td></tr><tr><td>Eccentricity at last stable orbit</td><td><i>e</i><sub>merge</sub></td><td>out of band</td></tr><tr><td>Signal to Noise Ratio</td><td>SNR</td><td>&lt; 50</td></tr></tbody></table>

## 2.6. Cosmic strings

Coordinators: Barry Wardell

Contributors: D. Chernoff, J. Wachter

**2.6.1. Description** Strings are effectively one-dimensional stress energy sources. If a network of strings is generated at early times then it can have many cosmological consequences including the production of gravitational waves (see [216] for a general review) that are potentially observable by LISA. String sources postulated by Grand Unified Theories (GUTs) have been ruled out by observations of the CMB [217, 218]. However, String Theory suggests new sources —fundamental strings and D-branes wrapped on small dimensions— in the context of certain inflationary string theory scenarios [219–232]. Macroscopic cosmic strings (CSs) are created and stretched to superhorizon scales by inflation (see [233, 234] for a review). Extensions to the Standard Model that introduce new symmetry breakings in-between the GUT and electro-weak scale can also produce viable strings [235, 236].

Irrespective of the detailed microscopic origin, these CSs may evolve to generate a network with some common features: after the Universe enters its radiation dominated phase the long, horizon-crossing strings begin to collide, break, reconnect and form small, sub-horizon scale loops. All string elements are dynamical, radiating gravitational waves and possibly other quanta [237–245]. Isolated loops, for example, radiate their entire rest mass energy and eventually disappear. Together the string elements generate a stochastic GW background [246–268] and GW bursts [252, 254, 258, 259, 264, 265, 268–274]. Detection and measurement of the string-generated gravitational waves by LISA will be informative for cosmology and high energy physics (see the whitepaper from the
LISA Cosmology Working Group for a summary [275]).

The physics of CSs is sensitive to (1) the set of fields to which the string couples, (2) whether the strings are global or local, (3) the ratio of the characteristic string width to curvature scale. Here, we assume that the strings are minimally coupled (they interact with each other but only radiate gravitational waves), local and well-described in the classical limit by the Nambu-Goto action. We refer to this as the “minimally coupled-string network”. There are additional possibilities but this defines a wide, interesting arena for this whitepaper.

Average properties of the minimally coupled string network are encapsulated in the Velocity One Scale (VOS) model [260, 276–278]. It turns out that all quantitative features depend primarily on the string tension $\mu_S$, or $G\mu_S/c^2$ in dimensionless terms. In particular, the total density in string components is parametrically small when the tension is small.

The VOS model is a valuable guide for forecasting the observations LISA may make. For specific numerical estimates below we assume minimally coupled local strings and adopt the following secondary parameters: intercommutation probability $\sim 1$ (intersecting field theory strings break and reconnect with probability of order unity [279] whereas string theory strings do so with smaller probability [280]), number of string species 1 (multiple species exist in realistic string constructions, for a review see [233]), fraction $\sim 0.2$ of long length strings chopped into loops of size $\sim 0.1$ of the horizon (for review of the small and large components inferred from simulation see [275, 281]) and rate of gravitational energy loss $dE/dt = \Gamma G\mu_S^2 c$ implied by dimensionless parameter $\Gamma = 50$ [247, 263, 282, 283]. Broadly speaking, changes to these adopted secondary values do not qualitatively change the network properties predicted by the VOS model. Many network properties are not included in VOS and have not yet been addressed in simulations. For example, isolated loops should evolve under the force of radiative backreaction (see [284–287]) but that process is not included in current numerical simulations and so it is difficult to accurately incorporate into statistical descriptions. These model-dependent, as opposed to parameter-dependent, uncertainties are important systematic deficiencies in our understanding and hard to quantify.

### 2.6.2. Tension limits

The most well-studied modern scenario involves Type IIB string theory and low tension strings produced at the end of brane inflation. Strings — not monopoles nor domain walls — are produced when a brane and anti-brane pair annihilate and initiate the Big Bang cosmology. The primary parameter, string tension, cannot currently be calculated a priori from theory. Instead we must turn to observations.

Empirical upper bounds on $G\mu_S/c^2$ have been derived from null results for experiments involving lensing [288–295], GW background and bursts [247, 252, 269, 270, 273, 296–305], pulsar timing [251, 253, 263, 266, 306–308], cosmic microwave background radiation
[217, 218, 309–319]. It has long been recognized [233] that all such bounds are model-dependent and typically involve observational and astrophysical uncertainties. Constraints from the CMB power spectrum rely on well-established gross properties of large-scale string networks and are relatively secure. Limits from optical lensing in fields of background galaxies rely on the theoretically well-understood deficit angle geometry of a string in spacetime but require a precise understanding of optical selection effects. Bounds from big bang nucleosynthesis rely on changes to the expansion rate from extra gravitational energy density but only constrain the strings formed prior to that epoch. Roughly speaking, these limits imply $G\mu_S/c^2 \lesssim 3 \times 10^{-8} - 3 \times 10^{-7}$ (see [252] for comparisons of limits). More stringent bounds on tension generally invoke additional assumptions [320]. Gravitational wave experiments [252, 274, 303, 321] can monitor the occurrence of bursts. In particular, the LVK set a bound $G\mu_S/c^2 \lesssim 4 \times 10^{-15}$ based on non-detection of assumed cusp-like bursts [322]. Long-term pulsar timing searches for a stochastic background have set the bound of $G\mu_S/c^2 \lesssim 1.5 \times 10^{-11}$ [263, 265, 266, 308]. In the future LISA may achieve limits as low as $G\mu_S \sim 10^{-17}$ for Nambu-Goto strings [275, 323].

### 2.6.3. Loop sources for LISA
The VOS model for the minimally coupled string network generates a loop size distribution weighted towards small sizes [264]. For string tensions $G\mu_S/c^2 \ll 10^{-7}$ the string loops are the most important elements of the network for GW science [275]. Small tensions imply weak GW damping. The undamped string is a non-linear oscillator with a fundamental period $T = \ell/(2c)$ and frequency $f = 1/T$. The Fourier transform of its motion yields power in all harmonics $nf$ for $n \ge 1$. A survey of the loop dynamics reveals large scale motions and distinctive small scale feature: cusps (infinitesimal bits of the string that move at the speed of light twice per fundamental period) and kinks (discontinuous changes of slope that perpetually circumnavigate the loop). Gravitational wave emission is sourced not only by the large scale oscillations but also by cusps and kinks. All long-lived loops are expected to possess cusps or kinks else they intercommute and produce kinks. Cusp-generated power decays with harmonic $n$ asymptotically $\propto n^{-4/3}$, kink-generated power $\propto n^{-5/3}$ and kink-kink collisions $\propto n^{-2}$ [267, 270, 272, 277]. At high frequencies, cusps dominate if they are present. Conversely, the period-averaged area of the sky illuminated by gravitational radiation increases from cusps to kinks to large scale modes.

There are three scenarios for LISA detections. (1) Loop decay creates a stochastic GW background from a large number of unresolved sources [275]. (2) Specific cusp or kink containing sources produce bursts of emission that stand above the general background [324]. (3) A few nearby loops, possibly associated with the Galaxy, produce emission that is strong, smooth and always on [281, 325–327].

### 2.6.4. Science with Cosmic Strings
(1) Interesting fundamental results are destined to emerge whether or not LISA detects evidence of gravitational radiation from strings. A
positive result for the stochastic background will allow the inference of the string tension;
a negative result will provide upper limits on the tension of any string component
that might be present [275]. If the network has a String Theory origin then either
determination helps guide progress towards a realistic model scenario for String Theory
that incorporates the Standard Model.

(2) A positive detection (either background, burst or nearby loop) is also fundamentally significant for cosmology because loops of macroscopic size are created during an epoch of inflation, supporting the inflationary paradigm [328–330]. The universe’s precise inflationary scenario remains a profound problem for cosmology and for fundamental physics. The almost scale-invariant density perturbation spectrum predicted by inflation is strongly supported by cosmological observations, in particular the cosmic microwave background radiation [218].

A negative detection is also very informative. The production of string-like structures
is a rather generic theoretical prediction whenever inflation does occur [236, 331]. A
negative result might be explained if the strings are unstable and/or couple to additional
fields that promote their decay. This sort of result will guide the search for models that
allow such interactions.
## 2.7. Beyond GR and beyond Standard Model sources

Coordinator: Paolo Pani, Helvi Witek

Contributors: N. Afshordi, R. Benkel, G. Bozzola, R. Brito, A. Cárdenas-Avendaño, E. Maggio, M. Okounkova, V. Paschalidis, C. Sopuerta

### 2.7.1. Introduction
BHs and compact binaries, and their GW emission have tremendous potential to probe for new physics beyond the Standard Model in the strong-field, nonlinear regime of gravity. LISA is likely to detect loud sources, such as MBHBs and EMRIs, which will allow us to test the nature of BHs, the validity of GR in the strong-field, highly-dynamical regime of gravity, or the presence of additional fundamental fields with unprecedented precision. These observations can help us address fundamental questions such as [82, 99, 102, 332–334]: (i) What is the nature of dark matter, and how can BH detections with LISA aid the search for new particles? (ii) What is the nature of gravity? Are there new fundamental fields and GW polarizations, as predicted by some extensions of GR and of the Standard Model? (iii) How do gravitational waves propagate over cosmological distances? (iv) Are the massive objects observed at galactic centers consistent with the rotating BHs predicted by GR? (v) Do exotic compact objects other than BHs and neutron stars exist in the universe? This complex topic in the context of LISA science is discussed in depth in the White paper of LISA’s Fundamental Physics Working Group [323]. Here we briefly summarize the most relevant sources for probing fundamental physics with LISA as a guide for the modelling of gravitational waveforms in GR and beyond.

### 2.7.2. Black holes and ultralight fields

#### 2.7.2.1. Description
Ultralight bosonic fields such as axions, axion-like particles or dark photons are predicted in several particle and theoretical physics models. A remarkable example is the string axiverse scenario, which predicts a multitude of axion-like particles emerging naturally from string theory compactifications [335]. These ultralight particles play a crucial role in diverse areas of physics and have been proposed (i) as a solution to the strong CP problem in quantum chromodynamics (QCD) [336], (ii) as compelling dark matter candidates [337, 338] and (iii) in cosmology [339]. Excitingly, we can employ BHs to search for (or constrain) ultralight bosons in a mass range that is complementary to traditional particle colliders or direct detection experiments [332, 340–344]. This surprising connection between BHs and particle physics is provided by the superradiant instability of BHs [340, 345–350]: low-frequency bosonic fields scatter off a rotating BH superradiantly, thereby extracting mass and angular momentum from the BH. Fields of mass-energy $\mu_B$ are efficiently confined in the vicinity of a BH with mass $M$ if the gravitational coupling $M\mu_B \lesssim 0.4$ [351–355], corresponding to $\mu_B \lesssim 10^{-16}(M/10^6 M_\odot)\text{eV}$. In this case they efficiently form a
bosonic condensate (“cloud”) around the BH. An alternative formation scenario involves accretion of such ultra-light fields onto BHs [356, 357] in the same mass range. If the bosonic field is complex, this process gives rise to hairy BHs [358–361].

The details of the cloud’s formation depends on the initial parameters such as the BH spin and the gravitational coupling between BH and bosonic field [362–364]. Once formed, the condensate dissipates by emitting a quasi-monochromatic GW signal [342–344, 362, 364–371]. The presence of boson clouds can also significantly affect the dynamics of binaries, e.g., through “dragging” of the cloud [372], dynamical friction [373] or tidal effects [374]. It also increases the eccentricity of comparable mass BH binaries and yields a GW phase shift [375]; In EMRIs the GW signal is modified relative to the vacuum case [376–379], and the presence of a secondary BH yields resonances [380–383]. The detection of EMRIs can be used to infer the bosons’ mass [384, 385].

2.7.2.2. *Expected source parameters* All “traditional” GW sources in the LISA band, including compact, comparable-mass BH binaries, EMRIs, and isolated spinning massive BHs, are potentially affected by bosonic clouds. Therefore, they are also good sources to act as cosmic laboratories for ultralight fields. Because the underlying mechanism only relies on the gravitational coupling, but is independent from the coupling of the bosons to the Standard Model of particle physics, they probe for all types of bosons, i.e., (pseudo-) scalars such as axion-like particles [335, 336], ultralight dark matter [337, 338], ultralight vector [386, 387] and tensor [371, 388] fields. LISA sources are particularly well suited for detecting or constraining ultra-light bosons in the mass range $\mu_B \in [10^{-19}, 10^{-15}]$ eV [343, 344, 371, 389, 390] (and even wider for massive spin-2 fields [388]), and they are suited for multi-wavelength searches in combination with ground-based instruments [391].

2.7.3. *Binary black holes as probes of the nature of gravity*

2.7.3.1. *Description* Does GR, our Standard Model of gravity, truly describe gravitational phenomena at all scales? It is expected to break down at high-energy scales as signaled by the presence of singularities inside BHs or at the Big Bang. At these scales a more complete theory of quantum gravity is needed that consistently combines gravity and quantum mechanics. However, GR cannot be quantized with standard approaches and it is not renormalizable. Therefore, GR (or a quantized version thereof) is not a viable candidate for quantum gravity. While a complete theory of quantum gravity remains elusive, most candidates predict similar extensions to GR such as higher curvature corrections or (non-minimal) coupling to new fields. BHs provide an ideal probe to search for such beyond-GR theories because, e.g., the presence of additional fields may endow BHs with scalar hair (or “charges”), thus violating the no-hair theorems of GR [392, 393]. More specifically, new fundamental fields arise in
several low-energy effective field theories of gravity [334], e.g., in the low-energy limit of quantum gravity, in the Horndeski class of scalar-tensor theories, tensor-vector-scalar theories [394], and in theories with quadratic curvature corrections. Such fields also arise in extensions of the Standard Model of particle physics, e.g., electric fields in mini-charged dark matter [395], primordial magnetic monopoles [396], darkly charged dark-matter [397], and the aforementioned bosonic clouds formed around BHs due to the superradiant instability.

In comparable-mass compact binaries, the coalescence of such “hairy” BHs generates additional scalar radiation that accelerates the inspiral and yields a GW phase shift. Furthermore, new polarization channels can exist in modified gravity theories. The detection or absence of such extra polarizations will be an important probe for new physics. In addition to modifications to the background solutions and the GW emission, modified theories of gravity may also change the physical properties of the gravitational waves once they are emitted, e.g., by changing the dispersion relation, the polarization and the way they interact with matter and with the detector [93, 333]. As modifications to the propagation of gravitational waves accumulate with the distance traveled, and the capability to put constraints on the mass of the graviton (Compton GW wavelength) scales with the chirp mass, comparable-mass binaries are the mose effective systems for measuring these effects [99].

EMRIs will provide an excellent probe of the multipolar structure of its primary object and, thus, test if the primary is consistent with the Kerr BH metric predicted by GR. GW measurements will reveal details of both the conservative (time-symmetric) and the dissipative (time-asymmetric) sectors of the gravitational theory. Additional degrees of freedom, such as dynamical scalar or vector fields, will introduce modifications to the motion of bodies, and additional sources of GW energy and angular momentum emission. Given that in most modified theories the gravitational field is described by a spin-2 metric tensor field and by additional fields [334], the interaction between matter and the new fields may give rise generally to an effective “fifth force”, leading to deviations from the universality of free fall [398], or in other words, to violations of the “strong” equivalence principle. How well a beyond-GR theory may be constrained depends on the relative PN order at which the correction enters and on the dimensions of the extra couplings [103, 399, 400]. In particular, theories with higher-order curvature invariants, such as the Gauss-Bonnet invariant [401–403] or the Pontryagin density [404] can be better tested with IMRIs and EMRIs [126, 134, 405].

2.7.3.2. *Expected source parameters* Both MBHBs and EMRIs can be employed to test gravity and the nature of compact objects with LISA [323]. For example, with nearly equal-mass BHs, one can perform null tests with inspiral-merger-ringdown consistency and BH spectroscopy, as well as searching for specific deviations with parametrized inspiral and ringdown waveforms. The beyond-GR modifications in the modelling of these signals (including PN theory, BH perturbation theory, effective-one-body
approaches, and numerical relativity) are discussed in Sec. 6.

With EMRIs one can constrain the multipolar structure of the primary object's spacetime with exquisite precision, thus testing whether the primary object is consistent with the Kerr BH metric predicted in GR. The parametrization of the waveforms will depend on the type of test carried out. For instance, when testing the geometry of the dark objects inhabiting galactic centers assuming GR, the parameters would describe the deviations from Kerr, e.g., multipole moments, tidal parameters or post-Kerr parameters. On the other hand, when testing GR, the parameters would describe the modifications of GR, e.g., additional coupling constants, extra dimensions length scale or higher-order corrections.

With both types of sources one can search for novel radiation channels due to extra polarizations or additional charges present in beyond-GR theories.

### 2.7.4. Testing the Kerr-hypothesis: BH mimickers and echoes

#### 2.7.4.1. Description
Exotic compact objects (ECOs) are horizonless objects which are predicted in some quantum gravity extensions of GR [406–410] and in the presence of exotic matter fields in the context of GR [102, 411–413]. The theoretical motivations for ECOs are the regularity of their inner structure and the overcoming of semi-classical puzzles such as that of the information loss [414–416]. These ideas have inspired a plethora of models including gravastars [417, 418], boson stars [419–422], wormholes [423–425], fuzzballs [416, 426] and others [427–432].

ECOs are classified in terms of their compactness, reflectivity and possible extra degrees of freedom related to additional fields [102, 433]. Two important categories are [434]: *ultracompact objects*, whose exterior spacetime has a photon sphere, and *clean-photon-sphere objects* (ClePhOs), so compact that the round-trip time of the light between the photon sphere and the object's surface is longer than the instability timescale of photon orbits. If the remnant of a merger is an ultracompact object, the ringdown signal differs from the BH ringdown at early stages and is dominated by the modified QNMs of the object [435, 436]. Conversely if the remnant of the merger is a ClePhO, the prompt ringdown is nearly indistinguishable from that of a BH because it is excited at the light ring [437]. The details of the object's interior appear at late times in the form of a modulated train of GW echoes [437–442]. The time delay between echoes depends on the compactness of the object and/or the energy scale of new physics [443], whereas the amplitude is related to the reflectivity of the object [102].

Many of these ECOs could be ruled out based on theoretical grounds. Horizonless compact objects are affected by an ergoregion instability when spinning sufficiently fast [444–447]. The endpoint of the instability could be a slowly spinning ECO [340, 448] or dissipation within the object could lead to a stable remnant [449, 450]. Furthermore, ultracompact horizonless objects might be generically affected by a light-ring instability
at the nonlinear level [448, 451, 452]. Current and future GW detectors will constrain
models of ECOs in almost all the regions of their parameter space [436, 453–455]. In
particular, searching for echoes in the post-merger signal of MBHBs with LISA [456]
will provide a clean smoking gun of deviations from the standard, “vacuum”, BH
prediction. Furthermore, EMRIs could constrain the reflectivity of the primary object
to unprecedented levels [457].

2.7.5. **Expected source parameters** The source parameters depend on the specific signal used to search for and constrain ECOs. The echo signal depends on the parameters of the remnant (in particular mass and spin), on its compactness, and especially on its effective reflectivity, which is zero for a classical BH [102, 458]. The reflectiviy can be generically a complex function of the frequency and of other remnant parameters. In an inspiral, besides the standard binary parameters, ECOs are characterized by anomalous multipole moments and nonvanishing tidal deformability, sharing in this case properties similar to those of BHs in modified gravity theories. In an EMRI, besides the different multipolar structure [459–470] and tidal deformability [471, 472] of the central object, a key parameter is again the (frequency-dependent) reflectivity [457].
### 3. Modelling requirements from data analysis

Ideally, waveform models should be infinitely accurate, evaluate instantly and be available in any format desired. In practice, none of these are achievable. The practical accuracy, efficiency and format requirements are set by the way LISA data is analyzed. This section starts with a brief overview of how LISA data analysis is expected to work. (For a more detailed description, see the LISA Data analysis whitepaper [473].) The remaining sections discuss how data analysis sets requirements on the accuracy, efficiency and formats of waveform models, providing the necessary framing for the discussion of waveform models in the rest of this whitepaper.

#### 3.1. Data analysis for LISA

Contributor: Tyson Littenberg

The foundation of gravitational wave data analysis is rooted in the conceptual simplicity of the measurement: observing relative changes in the separation between a collection of “proof masses” in free fall due to leading-order perturbations in the underlying spacetime metric, which propagate (effectively) uninhibited through the Universe. This is in stark contrast to, for example, electromagnetic observations, where the photons’ propagation from source to detector is influenced by intervening material (e.g., dispersion, scattering, absorption, reprocessing), and then undergoes complicated interactions with the instruments themselves (e.g., focusing optics, filters, diffraction, absorption by the detector) before registering as a signal. This is not to take away from the heroic effort and ingenuity required to develop the measurement system that is sensitive to the unfathomably small space-time perturbations themselves. However, given a detector that can achieve the necessary sensitivities, it is a tractable task to derive its response to incident gravitational waves from first principles.

As a result, GW data analysis methods have primarily developed around a *forward* problem, where the detector response is predicted from a hypothetical source and that prediction is then tested against the data [474]. There are notable exceptions, particularly in some searches for unmodelled GW transients in ground-based interferometer data [475], which approach the analysis as an *inverse* problem, starting with the observed data and working backward to solve for the input signal.

Any analysis is only as good as the models that go into it. The phenomenal sensitivity, accuracy and precision of GW observations is not achievable without highly accurate, coherent models for the gravitational waveforms themselves (the focus of this whitepaper), as well as the detector response and noise characteristics.

Under the assumption that the noise is Gaussian, the likelihood that the hypothetical
model, parameterized by $\theta$, would produce the observed data $\mathbf{d}$ is

$$p(\mathbf{d}|\theta) = \frac{1}{\sqrt{\det(2\pi\mathbf{C})}} e^{-\frac{1}{2}\mathbf{r}^\dagger\mathbf{C}^{-1}\mathbf{r}}, \quad (4)$$

where $\mathbf{r} = \mathbf{d} - \sum_i \mathbf{h}_i$ is a vector of all residual data samples after the discrete GW signals $\mathbf{h}_i$ have been subtracted, $\mathbf{C}$ is the noise covariance matrix $C_{ij} \propto \langle n_i | n_j \rangle$, and $\mathbf{n}$ is noise such that, in the absence of any GW signals, $\mathbf{d} = \mathbf{n}$. The likelihood is testing whether the residual is consistent with the ansatz that, in the absence of discrete gravitational waves, the data is Gaussian characterized by $\mathbf{C}$. Note the emphasis on discrete gravitational waves – a stochastic background of GW signals appears in the data model as a “noise” term, which modifies $\mathbf{C}$ with particular covariances that set it apart from instrument noise. See [476] for a detailed look at how the likelihood is derived, simplifications that produce slightly different forms of the equation in the literature, and a unified treatment of discrete and stochastic signals.

Also note that Eq. (4) does not prescribe a representation of the data. The game, as it were, is to represent the data in a way that minimizes the number of non-zero components of $\mathbf{C}^{-1}$ and thereby minimize the computational cost of evaluating the likelihood. It is the case under the assumption of stationary noise (i.e. $\mathbf{C}$ is constant over the observation time) that the discrete Fourier transform diagonalizes $\mathbf{C}$, which is why much of the GW analysis literature is based in the Fourier domain. For LISA, due to the long duration signals expected, the assumption of stationary noise will be dubious at best, and so analysis methods may trend towards other representations for the data (e.g., time-frequency methods with short Fourier transforms, discrete wavelet transforms, etc.), but the fundamental likelihood function remains the same.

The closest analog to the LISA analysis of individual sources is found in the analysis of ground-based interferometer data from the LVK collaboration, which are heavily reliant on waveform models. There are two major differences between LVK and LISA that limits the applicability of the analogy. First: At current detector sensitivities, the rate of detectable sources is such that they are still sparsely distributed through the data, with typical signals present in the most sensitive frequency band of the detectors for $O(10 \text{ s})$ at a rate of $O(\text{a few/week})$. LISA, on the other hand, will be signal-dominated, with tens of thousands of continuous Galactic sources overpowering the instrument noise below $O(3 \text{ mHz})$ and, depending on the rate and mass distribution of massive black hole mergers, several extremely high SNR mergers in band for $O(\text{weeks})$ to $O(\text{months})$, overwhelming any other contributions to the data stream. Thus, whereas the LVK searches are primarily “data mining” endeavors, sifting through a large volume of noisy data for rare and comparatively weak signals, LISA’s primary challenges are twofold: Source confusion due to the large number of sources simultaneously detectable; and model accuracy due to both the large number of waveform cycles over which models must stay phase coherent, and to contend with such high SNRs so as to not contaminate lower-amplitude sources with residual power.
The second key difference between the LVK collaboration and LISA experiences is the volume of the parameter space itself. Ground-based searches for compact mergers span a mass range for the components of O(1-100) M☉. While challenging, this mass range is small enough that precomputed grids of template waveforms covering the parameter space can be used when searching for candidates (see, e.g., [477]). The LISA parameter space for comparable-mass black hole mergers, for example, is several orders of magnitude larger, spanning O(10³−10⁸ M☉), both eliminating the possibility of using grid-based methods and expanding the range of possible mass ratios encountered by the analysis. To date, the most successful prototype LISA analyses have used stochastic sampling algorithms [478–480], still relying on template waveforms but using data-driven methods to concentrate waveform calculations in the high-probability regions of parameter space. While LVK analysis of compact mergers is hierarchical, with clear distinction between the “search” and “parameter estimation” steps, those two functions blur together for many prototype LISA pipelines. Stochastic sampling algorithms put more pressure on the computational efficiency of waveform calculations, since template generation is part of the analysis pipeline itself, as opposed to pre-computing and then reusing a (large) table of waveforms generated on a fixed grid.

Note that while 3rd generation ground-based GW detectors will trend towards higher event rates, signal durations and signal strengths, the LISA forecasted maximum signal strengths are uniquely in excess of signal-to-noise ratios ~ 10⁴ [481]. As a result, while the continued improvement of ground-based detectors and analysis methods naturally leads to evolution of the waveform models that directly benefit analysis of LISA data, the waveform development for 3rd generation detectors is necessary but not sufficient to fully achieve LISA’s potential. LISA puts unique pressure on the accuracy, breadth of parameter space and computational efficiency of waveform models.

It is also worth considering that the signal-processing part of the LISA science ground segment is divided between two paradigms: the so called “low latency” and “global” analyses. The global analysis is the joint fit to all GW signals in the data, and it is here where waveform accuracy is most important, in order to prevent mismatches with loud signals from contaminating weaker signals that are simultaneously present in the data stream. The global analysis is computationally intensive, and having efficient waveform generation tools is necessary for it to be tractable. However, the global-fit processing speed is more forgiving than that of the low-latency analyses, since the goal is to produce thorough source catalogs with a relatively relaxed release schedule on a O(monthly) to O(yearly) cadence.

For the low-latency analysis the trade-off is inverted, as computational speed and localization information are prioritized over all else. The low-latency pipelines will run ~daily, and will likely use a subset of the data (e.g. by bootstrapping based on the most recently completed analysis). The primary goal is to provide actionable information for joint multimessenger observations of transients. There are other functions of the low-latency pipelines, for example providing source-subtracted residuals to the instrument
team for assessment of detector performance etc., but under the purview of waveform generation time is of the essence. One thing to therefore consider at the architectural level of the waveform generation software is the need for tools that can be responsive to different demands on the speed/accuracy spectrum for processing sources, depending on the primary goal of the analysis.

## 3.2. Accuracy requirements

Coordinator: Deborah Ferguson and Maarten van de Meent

Contributors: M. Haney, R. O'Shaughnessy

Inaccuracies in waveform models affect the analysis of LISA data in three main ways. First, if the modelled waveforms do not sufficiently resemble the signals produced by Nature, this can hamper our ability to detect and identify sources in the data. Second, errors in the model will introduce some level of bias in the estimation of the source parameters, and could potentially masquerade as beyond-GR effects. Finally, if particularly loud sources are not perfectly subtracted from the data stream, their residuals can contaminate the searches for other sources. Below, we discuss the impact of these effects and how they lead to accuracy requirements for LISA waveform models.

### 3.2.1. Detection and identification
The impact of modeling errors on detection rates are fairly well understood in the case of a matched filter search for a single source in the data [482–484]. Suppose we have some waveform model, $h_{\text{model}}(\vec{\lambda})$, depending on some set of parameters $\vec{\lambda}$, and we are looking for some true waveform $h_{\text{true}}$ in the data, then the fitting factor is defined by

$$ \mathcal{F} = \max_{\vec{\lambda}} \frac{\langle h_{\text{true}} | h_{\text{model}}(\vec{\lambda}) \rangle}{\langle h_{\text{true}} | h_{\text{true}} \rangle \langle h_{\text{model}}(\vec{\lambda}) | h_{\text{model}}(\vec{\lambda}) \rangle}, \quad (5) $$

where $\langle \cdot | \cdot \rangle$ is the usual noise weighted inner product,

$$ \langle h_1 | h_2 \rangle = 4 \Re \int \frac{\bar{h}_1(f) h_2(f)}{S_n(f)} df, \quad (6) $$

with $S_n(f)$ the noise power spectral density of the detector [6, 7]. The fitting factor measures the effective loss in SNR due to using an imperfect model. Consequently, if the used models have a fitting factor $\mathcal{F}$ for a particular source, the maximum range at which such a source can be detected is reduced by a factor $\mathcal{F}$. Whether this has any impact on the detection rate depends on the type of source. Some LISA MBHBs sources are so loud that even the earliest (and therefore furthest) such events would be easily detectable [33], in which case the fitting factor of the used model has little to no effect on the detection rate. For other, quieter sources (such as EMRIs, SOBHBs, or GBs) the detection rate is more range limited, and a poor fitting factor $\mathcal{F}$ could lead to a reduction of the number of detections by a factor $\mathcal{F}^3$ (assuming a uniform
distribution of the source through a spatially flat universe). An appropriate norm for
what degree of loss of sources is deemed admissible for achieving LISA's science goals
needs to be established. An additional consideration here is that increasing the number
of unresolved sources could adversely affect the searches for other sources such as any
cosmological GW backgrounds [485, 486].

The above is valid for an idealized case, where the search is conducted using a continuous
bank of templates. In practice, a search would use a discretized template bank, meaning
that the effective fitting factor is increased due to the template spacing. Moreover, as
noted in Section 3.1, a fully coherent search of the LISA data seems infeasible. Instead,
LISA pipelines will most likely employ semi-coherent [487, 488] or stochastic [489, 490]
search strategies. In a semi-coherent search, template and data are both partitioned
into short segments, and the template-data overlap for each segment is maximized over
several extrinsic parameters. Waveforms used in such a search need only to stay phase-
coherent over the shorter segments. For example, it was shown in [488] that templates
with > 97% overlap accuracy over ≳ 10⁶ s will still be sufficient to detect > 50% of
EMRIs detectable with fully coherent, precise templates.

3.2.2. *Parameter estimation* LISA parameter inference nominally involves a joint multi-source fit for all available sources in the data [491–494]. For the purposes of the discussion below, we approximate this process as independent parameter inference for individual sources, resolving the source from detector noise and the confusion noise of all other signals in the data. Within that context, the impact of waveform sytematics on parameter inference has been historically estimated using well-understood analytic techniques; see, e.g., [495, 496]. The most frequently used ingredients in analytic waveform standards are the match or faithfulness *M*, which is the overlap between a signal *h*<sub>true</sub> and template *h*<sub>template</sub> maximized over the coalescence time and phase of the template; the Fisher matrix ⟨∂<sub>*a*</sub>h(λ)|∂<sub>*b*</sub>h(λ)⟩; and the inner product ⟨∂<sub>*a*</sub>h|δh⟩ between derivatives of the waveform and residuals δh between two signal approximations. The simplest and most conservative waveform accuracy standards for parameter inference are expressed in terms of limits on the mismatch 1 − *M*, which nominally must be ≲ 1/SNR<sup>2</sup> to avoid introducing systematic bias comparable to the statistical error for a single source [482, 483, 497]. For a population of sources, in principle systematic biases could stack in hierarchical population inference, and the most conservative threshold would be 1/(N<sub>s</sub>SNR<sup>2</sup>) where N<sub>s</sub> is the typical number of sources [498, 499]. However, this threshold is likely to be much too conservative for many applications, as the impact of systematic errors may be much smaller than features in the source population; see, e.g., [1].

Because the impact of systematic biases depends strongly on the nature and scale of
the bias relative to astrophysical features, general conclusions about systematic bias
cannot be drawn, and (barring negligible systematic error) must be assessed for each
science goal individually by performing hierarchical population inference. For example,
most SOBHBs in quasicircular orbits are in wide, relatively slowly-evolving orbits in the LISA band. For 90% of LISA-relevant SOBHBs (in population models consistent with current LVK observations), 2PN-accurate waveforms are sufficiently faithful for parameter estimation without biases [62]. At the other extreme, MBHBs will have extremely high amplitudes, with minute statistical error [492–494, 500]. Nominally, the very conservative accuracy thresholds described above would suggest a mismatch error target of order $1/(N_sSNR^2) \approx O(1/(10*(1000)^2)) \approx 10^{-7}$, for statistical errors to be small compared to statistical errors for the recovered population of MBHBs [496]. This threshold may be needed for applications that require joint inference on all massive MBHBs (e.g., tests of general relativity), but can be dramatically relaxed for any astrophysical interpretation of the MBHB population. Further studies will be necessary to establish specific requirements on waveform accuracy for each context in which parameter estimation is critical.

### 3.2.3. Contamination
Contamination effects occur when the residual from one source cannot be fully removed and impedes the interpretation of other, subdominant sources in the data. For such strong sources, systematic errors may be larger than statistical errors [496]. Though measurement error may have little impact on the astrophysical interpretation of the strong sources, the potentially substantial residuals produced by inaccurate models of their gravitational waves will introduce artifacts which contaminate downstream data analysis [491–494, 501]. For instance, for numerical relativity waveforms, insufficiently resolved grids can cause significant residuals even for simulations with precisely the same parameters as the observed signal [2]. The impact of errors arising from such imperfect modeling of strong sources can only be properly assessed with full joint hierarchical inference, using realistic models and contamination targets. Because of the diversity of sources that could be contaminated by residuals left by imperfect models for MBHB sources, much remains to be done to comprehensively assess key questions: How do modelling residuals from loud sources impact the detection and identification of other sources? Which sources are most affected? And what are the implications for the required modelling accuracy?

## 3.3. Efficiency considerations

Coordinator: Mark Hannam and Jonathan Thompson

Contributors: A. Chua and M. Katz

We shall discuss here the efficiency requirements imposed by search and inference for the two main classes of strong-field source that will be observable by LISA: MBHs and EMRIs. Relatively weaker-field sources such as GBs and SOBHs are a lesser concern; modeling of their signals (to the accuracy required for data analysis) is significantly easier, and computational cost is dominated by the need to attain a sampling resolution that is adequate for the typical duration, bandwidth and abundance of each source type.
Thus the efficiency requirements for these sources may be addressed on the analysis end with suitable approximations, or at the modeling–analysis interface where techniques such as those in Section 5 may be applied.

**3.3.1. Extreme-mass-ratio inspirals** EMRIs are the only LISA sources that combine the issue of strong-field complexity with that of long-lived signals, and thus they pose great challenges for both waveform modeling and data analysis. The main problem in EMRI search is *information volume*: the space of LISA-observable EMRIs is gargantuan, requiring $10^{30}-10^{40}$ templates to cover in a naive grid-based search [487].‡ This is unfeasible regardless of waveform efficiency. One possible solution relies on semi-coherent filtering, which essentially returns less informative templates that cover larger regions in parameter space, allowing extensive searches to be performed in a viable amount of time. Assuming that computing performance in the 2030s lies around the $10^2$-teraflop regime (achievable at present with ~ $10^3$-node clusters), semi-coherent segments of $\lesssim 10^6$ s would enable the analysis to be completed over the mission lifetime [487].

In practice, EMRI search will probably use semi-coherent filtering within stochastic algorithms, rather than template banks. Thus the estimates in [487] are very conservative, as stochastic searches are generally far more efficient than grid-based searches at the same number of template evaluations, even with a large multiplier (e.g., repeated runs) to ensure proper coverage of the parameter space. Note, however, that a potentially large number of additional templates will be required in the assessment of detection significance for each candidate source, although it remains unclear whether such an assessment will eventually be performed in the search or inference stage.

The efficiency requirements for inference are more straightforward, in that the analysis algorithms will not differ too much from those used in LVK parameter estimation. An estimated $10^6-10^9$ templates will be required per source posterior distribution. If we want to produce a posterior in less than 10 days we need waveforms that can be produced in $\lesssim$ 1ms. The estimated number of template evaluations per posterior also assumes that the prior regions for posterior sampling can be sufficiently localized to begin with; this is certainly possible in principle, but may require multiple search stages beforehand in a hierarchical approach.

**3.3.2. Massive black hole sources** Generally, MBH sources are expected to be detectable for ~week to hours prior to the merger. Due to this detection time frame and large signal-to-noise ratio of MBH binaries, coherent analysis over the full waveform template will most likely be employed for both search and parameter estimation. Therefore, from a waveform production perspective, these two analyses are roughly the same even if search and parameter estimation employ different schemes for achieving their various goals [502]. One way a search algorithm may differ from a parameter

<footnote>‡ As a point of reference, banks of $10^5-10^8$ templates are used to detect LVK BBH mergers [215].</footnote>
estimation algorithm for MBHs, in terms of the waveform generation, is in the use of higher-order harmonics. Recent work [503, 504] has shown the importance of higher-order modes in parameter estimation for LISA. However, as detailed in [502], searching over the dominant harmonic may be acceptable for initially and roughly locating sources throughout the high-dimensional parameter space. With this said, this is a matter of adding or removing modes (that we assume are available in a model), not altering the fundamental waveform generation method.

Given this idea that the waveform generation is similar for both search and parameter estimation for MBHs, the following discusses how efficiency requirements relate to the stages of waveform creation. As discussed in Sec. 5.1, there are generally two main parts to waveform creation in the context of LISA: sparse, accurate waveform calculations and a scaling method to achieve a full waveform from the sparse information. MBH waveforms can be generated beginning with sparse calculations of the amplitude and phase for each harmonic mode with an accurate waveform model. Depending on the specific analysis type, these sparse calculations are then upsampled to the desired search or parameter estimation settings. Producing accurate and upsampled waveforms *directly* from the accurate waveform generator is unnecessary and time-consuming: various methods allow for the upsampling of sparse, smooth functions in an accurate and much more efficient manner. This construction leads to two different waveform generation efficiency requirements. The first is the overall waveform (including upsampling). With similar requirements to EMRIs given in the section above, we expect to collect at least ~10⁶ MCMC samples for a given source posterior distribution. To accomplish this in ~10 days, we need to generate waveforms at a rate of ~1/s on a single central processing unit (CPU) core.

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

## 3.4. *Interface and data format requirements*

Contributors: Tyson Littenberg

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

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

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

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

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

In this section we outline the status of the above approaches and discuss the required
development needed to reach the accuracy requirements outlined in Sec. 3.2. Concerns
about the speed of waveform generation are addressed in Sec. 5. The status and requires
for LISA of these approaches can be found in Sec. 4.1 for numerical relativity, Sec. 4.2
for weak-field post-Newtonian and post-Minkowskian expansions, Sec. 4.3 GSF, Sec. 4.5
for EOB, and Sec. 4.6 for Phenomenological models. These sections focus on modelling
within GR; see Sec. 6 for a discussion on modelling in alternate theories of gravity.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Region of applicability of different approximation techniques for non-spinning quasi-circular binary BH inspiral. The 1PA region is a prediction derived from fitting NR data to an expansion of the form (14). The shaded regions indicate ranges within which the cumulative orbital phase-error is less than $\pi/4$ and $\pi/16$ radians, respectively. Recent direct calculations of post-adiabatic (1PA) waveforms suggest the 1PA region shown is over-optimistic for $\nu > 0.1$, but they have borne out the prediction for $\nu < 0.1$ [505]. The gray lines show the location of binaries with 10 (resp. 100 and 1000) orbits left before they reach the ISCO. This figure is reproduced from Ref. [506].</figcaption></figure>

## 4.1. Numerical relativity

Coordinators: Mekhi Dhesi, Deborah Ferguson, and Deirdre Shoemaker

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

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

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

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

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

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

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

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

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

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

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

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

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

Given the complexity of NR codes, cross-validation is an important component of code verification and it is by no means trivial that results obtained with different code bases agree. They employ, for example, different theoretical formulations of Einstein's equations for the evolution and initial data and different methods for the
<table><caption>**Table 8.** List of numerical relativity codes. We indicate if a code is open-source, if it has been used to produce public gravitational waveform catalogs, the formulation of Einstein's equation used (GHG: generalized harmonic, BSSN: Baumgarte-Shapiro-Shibata-Nakamura, CCZ4 / Z4c variants of the Z4 formulation, GCFE: generalised conformal field equations), if a code implements general relativistic hydrodynamics, and if it is capable of simulating compact binaries beyond general relativity. An asterisk indicates codes that are either (partially) based on the open-source Einstein Toolkit or are co-funded by its grant. Note this table was created jointly for this paper and Ref. [520].</caption><thead><tr><th>Code</th><th>Open Source</th><th>Public catalog</th><th>Formulation</th><th>Hydro</th><th>Beyond GR</th></tr></thead><tbody><tr><td>AMSS-NCKU [552, 613-615]</td><td>Y</td><td>-</td><td>BSSN/Z4c</td><td>-</td><td>Y</td></tr><tr><td>BAM [582, 616-618]</td><td>-</td><td>[524, 619, 620]</td><td>BSSN/Z4c</td><td>Y</td><td>-</td></tr><tr><td>BAMPS [621-623]</td><td>-</td><td>-</td><td>GHG</td><td>Y</td><td>-</td></tr><tr><td>COFFEE [624, 625]</td><td>Y</td><td>-</td><td>GCFE</td><td>-</td><td>Y</td></tr><tr><td>Dendro-GR [626, 627]</td><td>Y</td><td>-</td><td>BSSN/CCZ4</td><td>-</td><td>Y</td></tr><tr><td>Einstein Toolkit [628, 629]</td><td>Y</td><td>-</td><td>BSSN/Z4c</td><td>Y</td><td>Y</td></tr><tr><td>*Canada [366, 367, 630]</td><td>Y</td><td>-</td><td>BSSN</td><td>-</td><td>Y</td></tr><tr><td>*IllinoisGRMHD [631]</td><td>Y</td><td>-</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>*LazEv [515, 632]</td><td>-</td><td>[633-636]</td><td>BSSN/CCZ4</td><td>-</td><td>-</td></tr><tr><td>*Lean [637, 638]</td><td>Partially</td><td>-</td><td>BSSN</td><td>-</td><td>Y</td></tr><tr><td>*MAYA [639]</td><td>-</td><td>[639]</td><td>BSSN</td><td>-</td><td>Y</td></tr><tr><td>*NRPy+ [640]</td><td>Y</td><td>-</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>*SphericalNR [641, 642]</td><td>-</td><td>-</td><td>spherical BSSN</td><td>Y</td><td>-</td></tr><tr><td>*Spritz [643, 644]</td><td>Y</td><td>-</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>*THC [645-647]</td><td>Y</td><td>[619]</td><td>BSSN/Z4c</td><td>Y</td><td>-</td></tr><tr><td>*WhiskyMHD [648]</td><td>-</td><td>[649]</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>ExaHyPE [650]</td><td>Y</td><td>-</td><td>CCZ4</td><td>Y</td><td>-</td></tr><tr><td>FIL [651]</td><td>-</td><td>-</td><td>BSSN/Z4c/CCZ4</td><td>Y</td><td>-</td></tr><tr><td>GR-Athena++ [652]</td><td>Y</td><td>-</td><td>Z4c</td><td>Y</td><td>-</td></tr><tr><td>GRChombo [653-655]</td><td>Y</td><td>-</td><td>BSSN/CCZ4</td><td>-</td><td>Y</td></tr><tr><td>HAD [656-658]</td><td>-</td><td>-</td><td>CCZ4</td><td>Y</td><td>Y</td></tr><tr><td>Illinois GRMHD [659, 660]</td><td>-</td><td>-</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>MANGA/NRPy+ [661]</td><td>Partially</td><td>-</td><td>BSSN</td><td>Y</td><td>-</td></tr><tr><td>BH@H/NRPy+ [640, 662]</td><td>-</td><td>-</td><td>BSSN</td><td>-</td><td>-</td></tr><tr><td>MHDuet [663-665]</td><td>Y</td><td>-</td><td>CCZ4</td><td>Y</td><td>Y</td></tr><tr><td>SACRA [666-670]</td><td>-</td><td>[671]</td><td>BSSN/Z4c</td><td>Y</td><td>Y</td></tr><tr><td>SACRA-SFS2D [672, 673]</td><td>-</td><td>-</td><td>BSSN/Z4c</td><td>Y</td><td>-</td></tr><tr><td>SpEC [523, 674]</td><td>-</td><td>[521, 523, 675]</td><td>GHG</td><td>Y</td><td>Y</td></tr><tr><td>SpECTRE [676, 677]</td><td>Y</td><td>-</td><td>GHG</td><td>Y</td><td>-</td></tr><tr><td>SPHINCS_BSSN [678]</td><td>-</td><td>-</td><td>BSSN</td><td>SPH</td><td>-</td></tr></tbody></table>
GW extraction. They also employ different numerical techniques (e.g, pseudo-spectral or high-order finite differences discretization) and numerical implementations, e.g., for the grid structure. First systematic code comparisons and studies on the integration of NR with data analysis and analytical approximations have been conducted in Refs. [585, 587, 588, 592, 679, 680]. For example, the NINJA project [585, 592] focused on building a NR data analysis framework and to develop first injection studies. The initial study presented a sample of 23 NR waveforms of BBHs with moderate mass ratios $q \lesssim 4$, considering only the dominant the $\ell = |m| = 2$ GW mode. The follow-up study increased to statistics to 60 NR waveforms. While no cross-validation between the NR waveforms were performed, NINJA was crucial to identify technical and conceptual issues, that a hybridization with PN or EOB may be needed and that more than only the dominant GW mode may be necessary. The Samurai project [588] conducted a detailed cross-validation of gravitational waveforms obtained with five different NR codes. It focused on one BBH system, namely an equal-mass, non-spinning, quasi-circular (eccentricity $\lesssim 0.0016$) binary completing about six orbits before their merger. Focusing on the dominant $\ell = |m| = 2$ multipole, it was found that the waveforms' amplitude and phase agree within numerical error. It was also found that these NR waveforms were indistinguishable for SNR $\le 14$ and would yield mismatches of $\lesssim 10^{-3}$ for binaries with $M \sim 60M_\odot$ (using the anticipated detector noise curves at the time). The Numerical-Relativity-Analytical-Relativity (NRAR) collaboration [587] made important strides towards constructing more accurate inspiral-merger-ringdown waveforms by combining analytical and NR computations, and by defining new accuracy standards that would be needed for GW parameter estimation. The collaboration pushed cross-validation of NR codes to new frontiers by considering quasi-circular BBHs (with initial eccentricity $\le 2 \cdot 10^{-3}$) with both unequal mass (up to $q \le 3$) and moderately spinning BHs (up to $\chi_i \sim \pm 0.6$) completing about 20 GW cycles before merger. The study also included a binary of non-spinning BHs with $q = 10$, the highest mass ratio for quasi-circular binaries at the time. For the first time, the seven involved NR codes performed 22 targeted simulations to meet the required accuracy standards. They employed the same analysis code to estimate the uncertainties due to the numerical resolution and waveform extractions. The simulations exhibited a relative amplitude error in the $\ell = |m| = 2$ multipole of $\lesssim 1\%$ and a cumulative phase error of $\lesssim 0.25$rad. Finally, Ref. [680] presented a comparison of targeted simulations for the first GW event, GW150914, concluding that the waveforms were sufficiently accurate to effectively analyse LIGO data for comparable events. A comparison between different time evolution formulations of Einstein's equations, namely BSSN and Z4c, was performed in Refs. [552, 681].

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

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

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

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

Most NR simulations to-date have focused on quasi-circular inspirals because it is expected that any eccentricity present during a binary’s formation will have been emitted through GW radiation when it close to merger [192]. However, there are several astrophysical scenarios that predict eccentric binaries [696–703]. The NR community is beginning to explore the eccentric parameter space [596, 636, 704–717] as well as hyperbolic encounters [718–721].
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Parameter space coverage of public waveform catalogs in the quasi-circular limit. Shown are mass-ratio *q*, projections of the BH spins onto the orbital angular momentum (χ₁/₂z), and magnitude of spin-components orthogonal to the orbital angular momentum (χ₁/₂⊥) for waveform catalogs of different NR groups. Data as of June 2023.</figcaption></figure>

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

In summary, the BBH parameter space is currently covered as follows:
* Non-spinning BBH mergers with mass ratio $q \le 18$ [521, 523, 586, 587, 633, 634, 639, 726] — see [692, 727] for binaries with mass ratios $q = 64, 100$, and 128.

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

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

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

### 4.1.4. Environmental Effects

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

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

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

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

Finally, for a detailed discussion of environmental effects and the challenges they could present for precision tests of GR, such as dark matter environments, we refer the interested reader to the Fundamental Physics Whitepaper [323] and references therein.
**4.1.5. Challenges in NR** Improvements in hardware and numerical techniques continue to speed up simulations, yet some corners of the binary parameter space continue to demand challenging simulations that are costly in terms of runtime and computational resources. Three of the most pressing challenges facing NR are i) producing waveforms that cover the anticipated parameter space of unequal mass, highly spinning, and eccentric binaries, ii) producing long-lived waveforms with sufficient numbers of GW cycles, and iii) producing these waveforms at standards of accuracy set by LISA's sensitivity.

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

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

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

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

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

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

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

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

First, the explicit evolution schemes used in NR codes are subject to the Courant-Friedrichs-Levy (CFL) condition, which limits the possible time step to the dynamic
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Illustration of how numerical errors in a template may leave a residual in high-SNR signals. The left panel represents a BBH at high SNR (detector mass $M_{det} = 5 \times 10^6 M_{\odot}$, aligned spin of $\chi = 0.2$ on the larger BH at distance $D_L = 30$ Gpc and at inclination $\iota = 15^{\circ}$). The black line is the characteristic strain $(2f\tilde{h}(f))$ obtained at high numerical resolution from all spherical harmonic modes. The blue dashed line is the corresponding quantity obtained from a low resolution simulation, with the light blue line denoting the difference to the high-resolution waveform. The solid red line is a high resolution template containing only the $(l, m) = (2, 2)$ mode, with the light red line showing the difference to the black line. If this system were searched for with the low-resolution (or (2,2)-mode only) waveforms, then the shaded lines would remain as residuals after subtraction. The SNR of the residual due to resolution in the left panel is 342, clearly indicating the need for higher accuracy NR simulations. The right panel represents a more frequently expected system, where the SNR is 200 and $\iota = 0^{\circ}$. Here, the residual due to resolution has an SNR of 2.4. Figure adapted from [2].</figcaption></figure>

time of the smaller BH, $\Delta t \propto m_2$. As the mass ratio increases, the number of time-steps per unit $M = m_1 + m_2$ of evolution time will therefore increase as

$$f_{\text{CFL}} \sim \frac{q}{q_0}. \qquad (7)$$

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

$$f_{\Omega_{\text{low}}} = \left(\frac{\nu}{\nu_0}\right)^{-1} \left(\frac{\Omega}{\Omega_0}\right)^{-8/3}. \qquad (8)$$

Thirdly, LISA has higher accuracy requirements than current NR simulations. We denote by $f_{\tau}(A)$ the increase in computational cost that is needed to reduce the truncation error $\tau$ by a factor $A > 1$, i.e. $\tau \rightarrow \tau/A$.
A fourth cause for increased computational cost arises from the desire of achieving a certain overall phase error, which is important for GW waveform modeling. For longer simulations, phase accuracy must be preserved over more cycles. From PN theory, the orbital phase to merger scales as $\Phi \propto \nu^{-1}(M\Omega_i)^{-5/3}$. To achieve a fixed overall phase-error (say, 1 radian), the relative phase error must decrease $\propto 1/\Phi$, yielding an increase in computational cost by a factor

$$f_{\Phi} \sim f_{\tau} \left( \left( \frac{\nu}{\nu_0} \right)^{-1} \left( \frac{\Omega}{\Omega_0} \right)^{-5/3} \right). \qquad (9)$$

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

$$C = f_{\text{CFL}} f_{\Omega_{\text{low}}} f_{\tau} f_{\Phi}. \qquad (10)$$

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

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

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

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

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

For a more general discussion on waveform accuracy, see § 3.2.

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

Coordinator: Laura Bernard and Chris Kavanagh

Contributors: A. Antonelli, G. Faye, J. Garcia-Bellido, M. Haney, F. Larrouturou, A. Le Tiec, O. Long

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

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

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

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

As discussed below, much work has recently been dedicated to (i) push this accuracy to
<table><caption>**Table 9.** State-of-the-art of known PN results for both the conservative and dissipative dynamics as well as for the gravitational flux. Contrary to the text, everything is stated as absolute order. For example, the 6PN absolute order for the non-spinning flux in the table corresponds to the 3.5PN relative order results as stated in the text. Hereditary effects are only known as non-local contributions. An instantaneous expression can be obtained by performing a low-eccentricity expansion. SO and SS refer to spin-orbit and spin-spin interactions respectively.</caption><thead><tr><th rowspan="3">PN order</th><th colspan="4">Dynamics</th><th colspan="4">Dissipative flux</th></tr><tr><th rowspan="2">non-spinning</th><th colspan="3">spinning</th><th rowspan="2">non-spinning</th><th colspan="3">spinning</th></tr><tr><th>SO</th><th>SS</th><th>higher spins</th><th>SO</th><th>SS</th><th>higher spins</th></tr></thead><tbody><tr><td>0</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1.5</td><td>-</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2</td><td>√</td><td>-</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2.5</td><td>√</td><td>√</td><td>-</td><td>-</td><td>√</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3</td><td>√</td><td>-</td><td>√</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>3.5</td><td>√</td><td>√</td><td>-</td><td>√ (S<sup>3</sup>)</td><td>√</td><td>-</td><td>-</td><td>-</td></tr><tr><td>4</td><td>√</td><td>-</td><td>√</td><td>√ (S<sup>4</sup>)</td><td>√</td><td>√</td><td>-</td><td>-</td></tr><tr><td>4.5</td><td>*</td><td>√</td><td>-</td><td>√ (S<sup>3</sup>)</td><td>√</td><td>-</td><td>√</td><td>-</td></tr><tr><td>5</td><td>*</td><td>-</td><td>√</td><td>√ (S<sup>4</sup>)</td><td>√</td><td>√</td><td>-</td><td>-</td></tr><tr><td>5.5</td><td>*</td><td>-</td><td>-</td><td>√ (S<sup>5</sup>)</td><td>√</td><td>√</td><td>√</td><td>-</td></tr><tr><td>6</td><td></td><td></td><td></td><td>√ (S<sup>6</sup>)</td><td>√</td><td>√</td><td>√</td><td>√ (S<sup>3</sup>)</td></tr><tr><td>6.5</td><td></td><td></td><td></td><td>*</td><td>√</td><td>√</td><td>√</td><td></td></tr><tr><td>7</td><td></td><td></td><td></td><td>*</td><td>√</td><td></td><td></td><td></td></tr></tbody></table>

* means that only a partial result is known at those orders. At 5 and 5.5PN order, the results were obtained by the combination of PN traditional techniques with scattering amplitudes and self-force.

★ means that the dynamics is known at all leading order in the spin.

4.5PN order, (ii) include the effects of the spins of the compact objects (see also Table 9), (iii) generalize the results from circular to eccentric orbits, and (iv) explore the overlap and synergies with other approximation methods (PM approximation, small-mass-ratio approximation, effective one-body model). The reader is referred to the review articles [767, 778–785] and to the texbooks [786, 787] for more information.

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

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

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

4.2.3. Status

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

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

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

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

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

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

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

Finally, there have been some works trying to map the BMS (Bondi-Metzner-Sachs) flux-balance laws to equivalent PN results in the harmonic gauge [835]. By transforming the metric from harmonic coordinates to radiative Newman-Unti (NU) coordinates, they were able to obtain the mass and angular momentum aspects and the Bondi shear as a function of the quadrupole moment [836]. In particular, they rederive the displacement memory effect (see also [564]) and provide expressions for all Bondi aspects relevant to the study of leading and subleading memory effects [837].
4.2.3.2. Equations of motion and waveforms with spin. In the last decade, the effect of the proper rotation of the components of binary systems have been included in the orbital equations of motion to 4PN order, extending previous works dating back to the 70’s [838]. In the post-Newtonian terminology, the physical spin is rescaled as $S/c$. The leading spin interactions, the so-called spin-orbit interaction (SO), couple the mass $m_1$ of a particle with the spin of the other, hence those contributions scale as $\frac{Gm_1}{c^2} \times \frac{S_2}{c} \sim \mathcal{O}(\frac{1}{c^3})$ and are regarded as being of 1.5PN order. Meanwhile the first effect due to spin-spin interaction, scaling as $\frac{G}{c^2} \times \frac{S_1}{c} \times \frac{S_2}{c} \sim \mathcal{O}(\frac{1}{c^4})$, is quadratic in spin (SS) and arises at 2PN order. We refer to the successive subleading PN contributions to a given spin interaction (e.g. SO or SS) as the next-to-leading (NL), next-to-next-to-leading (NNL) and so on. Moving beyond the 2PN two-body dynamics including SO interactions [839], it becomes increasingly difficult to resort to an explicit model, where the bodies are described as small balls of fluid [840]. The most efficient strategy consists in adopting an effective point of view and considering that the two objects are point particles endowed with a classical spin.

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

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

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

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

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

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

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

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

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

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

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

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

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

The 3PM radiative contribution, which cancels a divergence in the 3PM conservative part, has been obtained in Refs. [942–948]. The radiative effects in PM expansions have been further explored in [949–952] and the radiative contributions to the scattering problem were obtained in a PN expansion [953–959].

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

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

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

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

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

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

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

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

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

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

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

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

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

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