# Science with the space-based interferometer LISA. V: Extreme mass-ratio inspirals

Stanislav Babak, Jonathan Gair, Alberto Sesana, Enrico Barausse, Carlos F.

Sopuerta, Christopher P.L. Berry, Emanuele Berti, Pau Amaro-Seoane,

Antoine Petiteau, Antoine Klein

► To cite this version:

Stanislav Babak, Jonathan Gair, Alberto Sesana, Enrico Barausse, Carlos F. Sopuerta, et al. Science with the space-based interferometer LISA. V: Extreme mass-ratio inspirals. Physical Review D, 2017, 95 (10), pp.103012. 10.1103/PhysRevD.95.103012. hal-01554417

HAL Id: hal-01554417

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

Submitted on 7 Oct 2025

**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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Science with the space-based interferometer LISA. V. Extreme mass-ratio inspirals

Stanislav Babak,¹ Jonathan Gair,² Alberto Sesana,³ Enrico Barausse,⁴ Carlos F. Sopuerta,⁵ Christopher P. L. Berry,³ Emanuele Berti,⁶,⁷ Pau Amaro-Seoane,⁵,⁸ Antoine Petiteau,⁹ and Antoine Klein⁴

¹Max Planck Institut fuer Gravitationsphysik, Albert-Einstein-Institut Am Muehlenberg 1,
D-14476 Golm, Germany

²School of Mathematics, University of Edinburgh, The King's Buildings, Peter Guthrie Tait Road,
Edinburgh EH9 3FD, United Kingdom

³School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT,
United Kingdom

⁴Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Université Paris 6 & CNRS, UMR 7095,
98 bis Bd Arago, 75014 Paris, France

⁵Institut de Ciències de l'Espai (CSIC-IEEC), Campus UAB, Carrer de Can Magrans s/n,
08193 Cerdanyola del Vallès, Spain

⁶Department of Physics and Astronomy, The University of Mississippi, University, Mississippi 38677, USA

⁷CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa,
Avenida Rovisco Pais 1, 1049 Lisboa, Portugal

⁸Institute of Applied Mathematics, Academy of Mathematics and Systems Science, CAS, Beijing 100190,
China Kavli Institute for Astronomy and Astrophysics, Beijing 100871,
China Zentrum für Astronomie und Astrophysik, TU Berlin, Hardenbergstraße 36, 10623 Berlin, Germany

⁹APC, Université Paris Diderot, Observatoire de Paris, Sorbonne Paris Cité,
10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

(Received 29 March 2017; published 31 May 2017)

The space-based Laser Interferometer Space Antenna (LISA) will be able to observe the gravitational-wave signals from systems comprised of a massive black hole and a stellar-mass compact object. These systems are known as extreme-mass-ratio inspirals (EMRIs) and are expected to complete ~$10^4–10^5$ cycles in band, thus allowing exquisite measurements of their parameters. In this work, we attempt to quantify the astrophysical uncertainties affecting the predictions for the number of EMRIs detectable by LISA, and find that competing astrophysical assumptions produce a variance of about three orders of magnitude in the expected intrinsic EMRI rate. However, we find that irrespective of the astrophysical model, at least a few EMRIs per year should be detectable by the LISA mission, with up to a few thousands per year under the most optimistic astrophysical assumptions. We also investigate the precision with which LISA will be able to extract the parameters of these sources. We find that typical fractional statistical errors with which the intrinsic parameters (redshifted masses, massive black hole spin and orbital eccentricity) can be recovered are ~$10^{-6}–10^{-4}$. Luminosity distance (which is required to infer true masses) is inferred to about 10% precision and sky position is localized to a few square degrees, while tests of the multipolar structure of the Kerr metric can be performed to percent-level precision or better.

DOI: 10.1103/PhysRevD.95.103012

I. INTRODUCTION

Gravitational waves (GWs) provide a means of gathering precious information otherwise beyond the reach of traditional electromagnetic astronomy. In particular, GWs may illuminate our understanding of the properties of black holes (BHs). The terrestrial Advanced LIGO [1] has recently observed GW signals from coalescing stellar-mass binary BHs, with two clear detections [2,3] and a probable third candidate [4,5]. These observations allowed estimation of the source parameters with high accuracy [5–7], giving new insight into their astrophysical formation [5,8] and allowing tests of general relativity (GR) [5,9,10]. Many more stellar-mass BH binaries are expected to be detected by LIGO (and by other terrestrial detectors such as

Advanced Virgo [11] and KAGRA [12]) in the next few years [5,13].

In addition to stellar-mass BHs, there is believed to be a population of massive BHs (MBHs), with masses in the range $10^5–10^9 M_\odot$, each lurking at the center of a galaxy [14–18]. Correlations between the mass of the MBH and other characteristics of the surrounding stars, such as the velocity dispersion $\sigma$ of the spheroidal component of the host galaxy (see, e.g., [19]) suggest a link between evolution of the MBH and its host galaxy [20–22].

Surrounding MBHs out to distances of a few parsecs, are nuclear star clusters of millions of stars [23]. In these innermost galactic regions, the density of stars easily exceeds $10^6 M_\odot$ pc$^{-3}$, and relative stellar velocities range
between ~100–1000 km s⁻¹. Here, mutual gravitational deflections between stars play a crucial role in determining dynamics [24], and their tidal disruption may contribute to increasing the mass of the central MBH [25–28]. Unlike stars, compact objects (COs; including stellar-mass BHs, neutron stars and white dwarfs) can avoid tidal disruption and approach the central MBH, radiating a significant amount of energy in GWs at low frequencies.

One of the main experimental challenges for ground-based detectors is seismic noise, which limits their sensitivity at frequencies ≲10 Hz, making them insensitive to GWs from MBH systems. However, space-borne interferometric GW detectors, such as the Laser Interferometer Space Antenna (LISA) [29], are free from the seismic noise. The technology behind LISA, an ESA-led mission expected to be launched by 2034, has been recently tested by the LISA Pathfinder experiment with outstanding results [30]. Previous work has investigated the scientific potential of LISA-like detectors for (i) MBH mergers and astrophysics [31]; (ii) stochastic backgrounds [32,33]; (iii) cosmography [34]; (iv) tests of general relativity [35,36]; and (v) ringdown tests of the nature of BHs [36,37]. LISA will also usher in the era of multiband GW astronomy, with stellar-mass binary BHs being detectable by LISA years to days before they reach the sensitivity window of ground-based detectors [38]. This would provide information on the formation mechanism of BH binaries [39–41], improve the precision of parameter estimation (including sky location) [42], and yield better constraints on putative deviations from GR [35]. In this paper we will focus on the physics and astrophysics of extreme mass-ratio inspirals (EMRIs) [43], i.e. systems comprised of stellar-mass BHs or other comparable mass COs orbiting around a MBH with mass M ~ 10⁴–10⁷ M☉.

As a consequence of their extreme mass ratio these systems inspiral slowly, completing ~10⁴–10⁵ cycles in LISA’s sensitive frequency range [44,45]. Therefore EMRIs are ideal signals to construct detailed maps of the background spacetime of MBHs [46–50], precisely estimate source parameters [51–53], perform tests of GR [50,54], and possibly detect the presence of gas around the central MBH [55–60]. Measuring the properties of a population of EMRI signals could additionally give us information on the mass distribution of MBHs [61] and their host stellar environments [43].

We examine in detail the scientific potential of EMRI observations with LISA, focusing on event rates and on parameter-estimation precision. There have been previous studies computing EMRI rates [62–64], but the astrophysical model employed in those calculations was a combination of simple power laws, and no attempt was made to quantify the uncertainties in that model. EMRI parameter-estimation studies have also been carried out [51,52], but only for a small sample of representative cases and not for a full astrophysical population. In this study we address both

of these shortcomings. We compute event rates for several different astrophysical models that were selected to quantify the main observational uncertainties, and we compute estimates of the parameter-estimation precisions for all the events in each population. Our results are computed for the first time considering a 2.5 Gm LISA detector with six laser links, which was proposed as the new mission baseline in the response to the ESA call in January 2017 [29].

The plan of the paper is as follows. We begin in Sec. II by discussing the assumed design of the LISA detector. In Sec. III we describe our astrophysical EMRI model and the related uncertainties. Section IV describes our EMRI waveform models and the parameter estimation calculation. We summarize our main results in Sec. V, and conclude by presenting possible directions for future research.

## II. LISA SENSITIVITY

The LISA baseline went through several stages of redesign in the past five years. Following the 2011 NASA drop-out, the classic LISA design was initially descoped to fit within the budget of an L-class ESA mission, leading to the New Gravitational-wave Observatory (NGO) design [65]. This new baseline was eventually selected as strawman mission in support of The Gravitational Universe [66], the science theme adopted by ESA for its L3 slot, scheduled for launch in 2034. Following the selection, in 2014 a Gravitational Observatory Advisory Team (GOAT) was appointed by ESA to consider a number of feasible options and issue a recommendation for a new design. The study considered a family of designs, featuring different choices for the arm length L, laser power, telescope diameter, mission duration and low-frequency noise level (see [31] for details).

Following the GOAT recommendation, the LISA Consortium answered the ESA call for missions by proposing the baseline outlined in [29]. The detector features a constellation of three satellites separated by L = 2.5 Gm and connected by six laser links. The output power of each laser is 2 W and their light is collected by 30 cm telescopes. The sky-averaged detector sensitivity can be written in analytic form as

$$S_n(f) = \frac{204S_n^{\text{acc}}(f) + 2S_n^{\text{loc}} + S_n^{\text{sn}} + S_n^{\text{omn}}}{3L^2} \times \left[1 + \left(\frac{2Lf}{0.41c}\right)^2\right], \quad (1)$$

where L is the arm length, and the noise contributions $S_n^{\text{acc}}(f)$, $S_n^{\text{loc}}$, $S_n^{\text{sn}}$ and $S_n^{\text{omn}}$ are due to low-frequency acceleration, local interferometer noise, shot noise and other measurement noise, respectively. The acceleration noise has been fitted to the level successfully demonstrated by the LISA Pathfinder [30] as
$$S_n^{\text{acc}}(f) = \left\{ 9 \times 10^{-30} + 3.24 \times 10^{-28} \left[ \left( \frac{3 \times 10^{-5} \text{ Hz}}{f} \right)^{10} + \left( \frac{10^{-4} \text{ Hz}}{f} \right)^2 \right] \right\} \left( \frac{1 \text{ Hz}}{2\pi f} \right)^4 \text{m}^2 \text{Hz}^{-1}, \quad (2)$$

whereas other contributions are set to

$$\begin{align}
S_n^{\text{loc}} &= 2.89 \times 10^{-24} \text{ m}^2 \text{Hz}^{-1}, \nonumber \\
S_n^{\text{sn}} &= 7.92 \times 10^{-23} \text{ m}^2 \text{Hz}^{-1}, \tag{3} \\
S_n^{\text{omn}} &= 4.00 \times 10^{-24} \text{ m}^2 \text{Hz}^{-1}. \nonumber
\end{align}$$

Besides the instrumental noise of Eq. (1), we also include a galactic confusion noise component, modeled by the fit

$$S_{\text{gal}} = A_{\text{gal}} \left(\frac{f}{1 \text{ Hz}}\right)^{-7/3} \exp\left[-\left(\frac{f}{s_1}\right)^{\alpha}\right] \times \frac{1}{2} \left[1 + \tanh\left(-\frac{f - f_0}{s_2}\right)\right]. \quad (4)$$

The overall amplitude of the background $A_{\text{gal}} = 3.266 \times 10^{-44} \text{ Hz}^{-1}$ depends on the astrophysical model for the population of white dwarf binaries in the Galaxy. Here we have used the same model as in [29]. The power law $f^{-7/3}$ is what we expect from a population of almost monochromatic binaries. The exponential factor comes from removal of the loud signals standing above the confusion background, while the last term takes into account that all Galactic binaries can be resolved and removed above some frequency $f_0$. For the assumed two-year observation period, the fitting parameters appearing in the above expression for $S_{\text{gal}}$ have the values: $\alpha = 1.183$, $s_1 = 1.426 \text{ mHz}$, $f_0 = 2.412 \text{ mHz}$, $s_2 = 4.835 \text{ mHz}$.

The LISA design is most sensitive at millihertz frequencies, making it well-purposed for observing EMRs.

### III. ASTROPHYSICAL EMRI MODEL

The expected EMRI rate depends on several astrophysical ingredients:

(i) The MBH population in the accessible LISA mass range, $M \in [10^4, 10^7] M_{\odot}$, the redshift evolution of their mass function, and their spin distribution;

(ii) The fraction of MBHs hosted in dense stellar cusps, which are the nurseries for EMRI formation;

(iii) The EMRI rate per individual MBH, and the mass and eccentricity distribution of the inspiralling COs. In the following subsections we consider these ingredients in turn, presenting the astrophysically motivated prescriptions used in this work, before combining them in Section III D.

#### A. MBH population

We consider here two population models that are intended to bracket current uncertainties in the MBH mass

<figure><img src="image_2.png" /><figcaption>FIG. 1. MBH density mass function $dn/d\log_{10} M$ for the self-consistent model popIII at redshift 0 (solid), 1 (long dashed), 2 (short dashed) and 3 (dotted). The approximation provided by Eq. (5) is shown as a thin straight black line. Also shown in brown is the redshift-independent pessimistic mass function as given by Eq. (6). The shaded area represent constraints from Shankar et al. [72] (light orange) and Shankar [73] (green).</figcaption></figure>

function at the low mass end (cf. Fig. 1). The first one is Model popIII, as investigated in Klein et al. [31]. This is a self-consistent model for MBH formation and cosmic evolution developed in [67–70], and assumes light MBH seeds from population III (popIII) stars [71], while accounting for the delays between MBH and galaxy mergers. The model successfully reproduces several galaxy and MBH mass function properties, and it is consistent with observational constraints on the MBH mass function [72,73]. The predicted MBH mass function in the relevant range can be approximated as

$$\frac{dn}{d \log M} = 0.005 \left( \frac{M}{3 \times 10^6 M_{\odot}} \right)^{-0.3} \text{Mpc}^{-3}, \quad (5)$$

almost independent of redshift, as shown in Fig. 1. We label this mass function “Barausse12”.

Following Gair et al. [61], we also consider a more conservative model with a redshift-independent mass function of the form

$$\frac{dn}{d \log M} = 0.002 \left( \frac{M}{3 \times 10^6 M_{\odot}} \right)^{0.3} \text{Mpc}^{-3}. \quad (6)$$

In this case, the MBH mass function increases with mass at the low-mass end, and it is therefore less favorable for EMRI events falling in the LISA band. This is a purely
phenomenological model, which does not come from a self-consistent MBH evolutionary scenario, but is still consistent with current observational constraints on the MBH mass function. We label this mass function “Gair10.”

The EMRI rate and expected signal also depend on the spin parameter a of the central MBH. The popIII model self-consistently follows the spin evolution of MBHs through accretion and mergers. We find that most MBHs in the LISA range have near maximal spins. This is because, although MBH seeds start with a random spin distribution, they need to accrete enough mass to get into the LISA band. At small masses, the MBHs always align with the accreted material (the MBH angular momentum is always smaller than the disk angular momentum in our model [68]). The distribution starts to extend to lower a for higher masses, when the MBH spin becomes larger than the typical disk angular momentum [74]. As a result, the MBHs do not always align with the accreting material, and spindown is possible. However, this effect becomes appreciable only at $M \approx 10^7 M_\odot$ [68,74]. We assume a maximum dimensionless MBH spin parameter $a = 0.998$, with a median value around $a = 0.98$. Since most MBHs have high spins in our default model, we label it “a98.” For the sake of comparison, we also consider two alternative models; one with a flat spin distribution $0 < a < 1$, labeled “aflat,” and one with nonspinning MBHs, labeled “a0.”

## B. Stellar cusps surrounding MBHs

A necessary condition for EMRI formation is the presence of a cusplike distribution of stellar objects surrounding the MBH. It has generally been assumed that MBHs are immersed in a Bahcall-Wolf stellar cusp with density profile $\rho(r) \propto r^{-7/4}$, which is the steady state solution for a distribution of stars in the sphere of influence of a massive object [75]. However galaxies merge, and so do the MBHs they host. MBH binaries destroy stellar cusps, carving a low density core [69,70,76] which is unsuitable to the formation of EMRIs. One of the main advantages of using a semianalytic MBH evolution model is that we are able to track the MBH merger history implementing a simple prescription that takes into account in a self-consistent way cusp disruption following MBH binary mergers.

### 1. Cusp erosion and regrowth

To understand the impact of a merger we must estimate the time $t_{cusp}$ taken for a cusp to reform.

We assume that each MBH binary with mass $M = M_1 + M_2$ is embedded in an isothermal sphere, defined by a density profile [77]

$$ \rho(r) = \frac{\sigma^2}{2\pi G r^2}, \qquad (7) $$

where $\sigma$ is the one-dimensional velocity dispersion. We further assume that the MBH binary carves a core of

constant density and size $r_c$ in the center of the stellar system. The mass deficit due to a flat core of size $r_c$ is given by

$$ M_d = \frac{4 \sigma^2 r_c}{3 G}. \qquad (8) $$

Thus, $r_c$ can be estimated once $M_d$ is known. The mass deficit must equal the mass displaced by the MBH binary on its way to coalescence, and is estimated as [69,70]:

$$ M_d = 0.7 M q^{0.2} + 0.5 M \ln\left(\frac{r_h}{r_{GW}}\right) + 5 M \left(\frac{V_k}{V_{esc}}\right)^{1.75}. \qquad (9) $$

Here, $q = M_2/M_1 \le 1$ is the mass ratio of the MBH binary, $r_h$ is the binary hardening radius, $r_{GW}$ is the radius at which GW emission dominates over stellar hardening, $V_k$ is the GW kick and $V_{esc} \approx 5\sigma$ is the typical escape velocity from the stellar bulge [78]. To make use of Eq. (9), we need an estimate of $r_h/r_{GW}$. Here $r_h$ is the hardening radius, the separation at which the specific binding energy of the binary is equal to the average specific kinetic energy of the surrounding stars [79,80], given by

$$ r_h = \frac{GM_2}{4\sigma^2}, \qquad (10) $$

where $M_2$ is the secondary’s mass (see e.g. [81]). The distance $r_{GW}$ represents the separation at which the MBH binary evolution switches from being stellar hardening dominated to be GW driven. It can therefore be computed by finding where the three-body scattering hardening rate $(dr/dt)_*$ becomes equal to the GW shrinking rate $(dr/dt)_{GW}$. The latter given by the standard quadrupole formula [44], and the former can be written as [81]

$$ \left(\frac{dr}{dt}\right)_* = \frac{HG\rho_*}{\sigma} r^2. \qquad (11) $$

Here $H \approx 15$ is a dimensionless hardening rate and the stellar density $\rho_*$ is evaluated at the influence radius of the binary $r_i = GM/\sigma^2$ [82]. For the isothermal sphere this gives

$$ \rho_* = \frac{\sigma^6}{2\pi G^3 M^2}. \qquad (12) $$

Combining everything together and assuming circular binaries, one gets:

$$ \frac{r_h}{r_{GW}} \approx 0.178 \frac{c}{\sigma} \frac{q^{4/5}}{(1+q)^{3/5}}. \qquad (13) $$

For a given MBH mass, if we know that there was a merger with a given $q$, we can substitute Eq. (13) into Eq. (9), to obtain the mass deficit, and use this in Eq. (8) to solve for $r_c$ and obtain the extent of the core. Once $r_c$ is known, the relaxation time for an isothermal sphere is given by [77]
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2. Left panel: Cusp regrowth time $t_{\text{cusp}}$ as a function of the total MBH binary mass. Solid, long-dashed and short-dashed curves are for $q = 1, 0.1, 0.01$ respectively. Red curves assume $V_k = 0$ whereas blue curves assume $V_k/V_{\text{esc}} = 0.6$. Right panel: Mass deficit normalized to $M$ as a function of binary mass ratio for $M = 10^5 M_{\odot}$ (short dashed), $M = 10^6 M_{\odot}$ (long dashed), and $M = 10^7 M_{\odot}$ (solid). Blue and green dots are mass deficits computed by Khan et al. [84].</figcaption></figure>

$$t_{\text{relax}} = \frac{5}{\ln \Lambda} \left( \frac{\sigma}{10 \text{ km s}^{-1}} \right) \left( \frac{r_c}{1 \text{ pc}} \right)^2 \text{Gyr}, \quad (14)$$

where $\ln \Lambda \approx 10$ is the Coulomb logarithm [77]. The cusp regrowth time is then [83]

$$t_{\text{cusp}} = 0.25 t_{\text{relax}}. \quad (15)$$

This can be expressed in terms of $M$ and $q$ only if we specify an $M-\sigma$ relation to eliminate the $\sigma$ dependence. We use the best fit of Gültekin et al. [18] as our default model:

$$M = 1.53 \times 10^6 \left( \frac{\sigma}{70 \text{ km s}^{-1}} \right)^{4.24} M_{\odot}. \quad (16)$$

Combining Eqs. (8)–(15), it is possible to approximate $t_{\text{cusp}} \propto M^{1.29}$ if we ignore the $r_h/r_{\text{GW}}$ term in Eq. (9). The dependence on the mass ratio is mild. Results for $t_{\text{cusp}}$ are shown in the left panel of Fig. 2. The red curves are for $V_k = 0$, whereas the blue ones assume $V_k = 0.6V_{\text{esc}}$. If we ignore the $V_k$ contribution, we can fit the cusp regrowth time as

$$t_{\text{cusp}} \approx 6 \left( \frac{M}{10^6 M_{\odot}} \right)^{1.19} q^{0.35} \text{Gyr.} \quad (17)$$

The slightly weaker dependence on $M$ than the initial approximation is due to the $r_h/r_{\text{GW}}$ term in Eq. (9). Typical cusp regrowth time scales are a significant fraction of the Hubble time for equal-mass binaries with total mass $10^6 M_{\odot}$, whereas they tend to become unimportant for

lower mass MBHs (generally less than 1 Gyr for a $10^5 M_{\odot}$ MBH).

Further core scouring following significant kicks will make these time scales a factor of 2 longer. For typical kick velocities of few hundred km s⁻¹ we find that the EMRI rate drops by a factor of ~2 due to a combination of MBH ejections from low mass halos and prolongation of cusp regrowth time scales.

The adopted $M-\sigma$ relation has a significant impact on $t_{\text{cusp}}$. We will therefore explore different prescriptions. As mentioned above, our default model employs the $M-\sigma$ relation of Gültekin et al. [18] (labeled “Gultekin09”) which gives $t_{\text{cusp}} \approx 6$ Gyr for a Milky Way-like MBH. We also consider two alternatives; a pessimistic model from Kormendy and Ho [19] (labeled “KormendyHo13”) which gives $t_{\text{cusp}} \approx 10$ Gyr for a Milky Way-like MBH, and an optimistic one [85] (labeled “GrahamScott13”) which gives $t_{\text{cusp}} \approx 2$ Gyr for a Milky Way-like MBH.

For the sake of completeness, we also tried a model based on Shankar et al. [86], which claims that the observed $M-\sigma$ relations are fundamentally biased and that the intrinsic one has a lower normalization. We found this to make little difference in practice for EMRI rates, and do not present results based on this model.

To verify our simple model, we performed a series of sanity checks. First, for a Milky Way-like MBH, Eq. (14) implies $t_{\text{relax}} \approx 10^{11}$ yr, which is consistent with Fig. 1 of [87]. Moreover, $M_{\text{d}}$ given by Eq. (9) is consistent with the results of full N-body simulations by Khan et al. [84],
<figure><img src="image_5.png" /><figcaption>FIG. 3. Cusp regrowth effect for the popIII model. Left panel: The average differential number of mergers per unit redshift (i.e. Eq. (18) integrated over $q$) $dN_m/dz$ experienced by each individual MBH of mass $\log_{10} M = 4.5, 5, 5.5, 6, 6.5$ from darker-thicker to lighter-thinner. Center panel: The solid curves are the values of $N_m(M, z)$ given by Eq. (20), and the dashed curves are the corresponding probabilities of retaining a cusp given by Eq. (21). Right panel: The differential number of MBHs $dN/dz$ across the Universe in the three different mass bins that are potential EMRI hosts, either ignoring cusp disruption (solid lines) or taking it into account (dashed lines).</figcaption></figure>

as shown in the right panel of Fig. 2. Here, each of the red lines shows the mass predicted by our simple model as a function of $q$ (different line styles refer to different MBH masses). The blue and green dots are mass deficits computed by Khan *et al.* [84] at the end of their simulations. The blue dots are mass deficits within $1.5r_i$, where $r_i$ is the MBH binary influence radius, whereas the green dots are mass deficits within $3r_i$. The mass deficit in those simulation saturate between $2r_i$ and $3r_i$. The cores predicted by our simple model are $\approx 1.5r_i$. In the simulations, however, the MBH binaries do not evolve all the way through coalescence. Small $q$ binaries, in particular, are stopped at an earlier stage of the evolution, because the simulations are more time consuming. This is the likely explanation of the steeper mass-ratio dependence of the simulation results with respect to our models. Overall, the analytical mass deficits and the results of the simulations agree to within a factor of 2.

## 2. Fraction of MBHs hosted in stellar cusps

To compute the fraction of MBHs that reside in cusps versus those in cores, we need to convolve the MBH number density $dn/dM$ (ignoring the spin dependence for the moment) with the number density of mergers per unit mass, mass ratio and redshift $d^3n_m/dMdzdq$, and the cusp regrowth time $t_{cusp}(M, q)$ given by Eq. (17).

First, we assume that MBHs do not grow appreciably in mass in the redshift range of interest (mostly $z < 2$ for LISA). Although this might well be a crude approximation, it simplifies the model. From our semianalytic MBH evolution model [67–70], we extract the distribution $d^3n_m/dMdzdq$, which is the differential number density (per Mpc$^3$) of mergers with mass ratio $q$ undergone by a MBH of a given mass $M$ at redshift $z$. The quantity of interest is $p_0(M, z)$, the probability that a MBH of mass $M$

observed at redshift $z$ had zero mergers within its cusp regrowth time $t_{cusp}(M, q)$, given by Eq. (17). We can define the quantity $d^2N_m/dzdq$ as

$$ \frac{d^2N_m}{dzdq}(M, z, q) = \frac{d^3n_m}{dMdzdq} \left(\frac{dn}{dM}\right)^{-1}. \quad (18) $$

This is the (mean) differential merger rate for an individual MBH with mass $M$, i.e. the number of mergers a MBH of a given mass $M$ has undergone between redshift $z$ and $z + dz$ and with mass ratio in the range $q$ and $q + dq$. The integral over $q$ of Eq. (18) is represented in the left panel of Fig. 3; MBHs in the mass range of interest for LISA generally experience between 0.1 and 1 merger per unit redshift since $z = 6$. For each mass ratio, we can then define a critical redshift $z_{cusp}(M, q)$ by solving the implicit equation

$$ t_{cusp}(M, q) = \int_{z}^{z_{cusp}(M, q)} dz' \frac{dt}{dz'}, \quad (19) $$

where $t_{cusp}(M, q)$ is computed using Eq. (14) and Eq. (15); if a MBH had suffered a merger between $z$ and $z_{cusp}(M, q)$, there would be no cusp. We can then compute the mean number of mergers $N_m(M, z)$ experienced by an individual MBH of mass $M$ observed at redshift $z$ in its cusp regrowth time as

$$ N_m(M, z) = \int dq \int_{z}^{z_{cusp}(M, q)} dz' \frac{d^2N_m}{dz'dq}(M, z, q). \quad (20) $$

Assuming Poissonian statistics for the mergers, the probability that a MBH of mass $M$ and redshift $z$ did not suffer a merger within its cusp regrowth time is<fnref n="1" />

<footnote n="1">¹A Poissonian probability distribution is strictly speaking valid only for rare, statistically independent events with a constant rate per unit time. Nevertheless, one can easily show that this equation holds also for events with nonconstant rate.</footnote>
$$p_0(M, z) = \exp[-N_m(M, z)]. \quad (21)$$

We apply to each MBH a probability $p_0(M, z)$ of being hosted in a stellar cusp, therefore being a suitable candidate for capturing an EMRI. $N_m(M, z)$ is shown in the central panel of Fig. 3 for different MBH mass values. Despite the similar number of mergers across the mass spectrum (left panel), low mass MBHs observed at any $z$ are extremely unlikely to have undergone a merger within their short cusp regrowth time scale (see left panel of Fig. 2), and their probability of being hosted in a stellar cusp is of order unity. The opposite is true for massive MBHs which reside in galaxies with much longer cusp regrowth time scales, and are likely to be hosted in a low-density stellar core.

If the distribution of MBHs in the Universe is described by a mass function (now including spin) of the form $d^3N/(dMdzda)$, then we can construct an effective MBH mass function for MBHs which could be potential EMRI hosts:

$$\left(\frac{d^3N}{dMdzda}\right)_{\text{eff}} = \frac{d^3N}{dMdzda} p_0(M, z). \quad (22)$$

The right panel of Fig. 3 shows this distribution integrated in spin and in different mass bins. It is clear that the number of potential EMRI hosts is severely suppressed only for $M > 10^6 M_\odot$.

### C. EMRI rate per MBH and properties of the stellar-mass BH

Finally, we need to specify the rate $R_0$ at which COs are captured by the central MBH, and define the properties of their orbits. The CO capture rate by MBHs has been investigated extensively in the literature, taking into account the effect of mass segregation [83], resonant relaxation [88], relativistic corrections [89], central MBH spins [90] and initial density profiles of the COs [91].

Our starting point is the intrinsic rate from Amaro-Seoane and Preto [83], which accounts for the effect of mass segregation:

$$R_0 = 300 \left( \frac{M}{10^6 M_\odot} \right)^{-0.19} \text{Gyr}^{-1}. \quad (23)$$

This is useful scaling relation; however, it has been calibrated for Milky Way-like galaxies, and care must be taken when extrapolating to other systems. In particular, this rate was calculated assuming a steady-state stellar environment surrounding the (growing) MBH which often cannot be achieved, especially for low-mass MBHs. Moreover, Eq. (23) only describes the EMRI rate: it does not include direct plunges. COs can be scattered onto nearly radial orbits, directly plunging into the MBH without emitting a significant GW signal. Although such systems are lost as GW sources, they do contribute to the growth of the MBH. The ratio of plunges to EMRI

depends mostly on the steepness of the density profile of the CO population. Compared to EMRIs, plunges are typically scattered into the MBH from much greater distances, so that a flatter density profile results in a larger plunge-to-EMRI ratio. For example, Merritt [91] considered two different CO distributions around MBHs of $10^6 M_\odot$ and $4 \times 10^6 M_\odot$, and found that while the EMRI rate varied within a factor of 2, remaining consistent with Eq. (23), the number of plunges per EMRI, $N_p$, went from being less than one for the steeper density profile, to be more than 50 for the shallower one. A recent study including a single population of compact objects found more than 100 plunges per EMRI [92]. A proper computation of EMRI rates in an astrophysical context would require N-body simulations starting from realistic initial conditions, spanning a wide range of MBH masses and of their surrounding stellar distribution properties. This is a challenge that goes beyond the scope of this paper, and in the following we develop a simple model to quantify the impact of nonstationary CO feeding rates and direct plunges on the astrophysical EMRI rates.

The parameter $N_p$ introduced earlier can vary between zero and $\sim 10^2$. Using Eq. (23), the total mass accretion rate for the MBH is given by

$$\begin{aligned} \dot{M} &= m R_0 (1 + N_p) \\ &= 3000 (1 + N_p) \left( \frac{m}{10 M_\odot} \right) \left( \frac{M}{10^6 M_\odot} \right)^{-0.19} M_\odot \text{Gyr}^{-1}, \end{aligned} \quad (24)$$

where $m$ is the characteristic mass of the CO. There are two problems that arise from implementing Eq. (24), which are exacerbated for low MBH masses. Consider for example $M = 10^5 M_\odot$, $m = 10 M_\odot$ and $N_p = 10$. First, according to this prescription, the MBH would double its mass in only 2 Gyr, accreting more than five times its initial mass in a Hubble time. Therefore, accreting COs at the rate given by Eq. (23) would be inconsistent with the existence of $M = 10^5 M_\odot$ MBHs [93]. Second, such a high accretion rate implies an astrophysically implausible supply of COs to the MBH. Assuming a standard Salpeter mass function [94,95], only about 0.3% of stars have a mass $m_* > 30 M_\odot$. Assuming those end their life as COs of $m \approx 10 M_\odot$ [96,97], then we can estimate that about 3% of the total stellar bulge is indeed composed by COs. Within the sphere of influence of the MBH the enclosed mass in stars is $M_* = 2M$, and therefore the mass in remnant BHs is about $M_{CO} = 0.06M$. The CO content within the sphere of influence of the MBH would therefore be depleted in a time

$$\begin{aligned} t_d &= \frac{M_{CO}}{\dot{M}} = \frac{0.06M}{m R_0 (1 + N_p)} \\ &= \frac{20}{1 + N_p} \text{Gyr} \left( \frac{m}{10 M_\odot} \right)^{-1} \left( \frac{M}{10^6 M_\odot} \right)^{1.19}. \end{aligned} \quad (25)$$
This can be compared to the relaxation time defined by Eq. (14), where we can substitute $r_c$ with the influence radius $r_i \approx 2GM/\sigma^2$ of the central MBH. By using the $M-\sigma$ relation of Gültekin et al. [18] and assuming $\ln \Lambda = 10$, we can compute the ratio of the two time scales:

$$ \frac{t_d}{t_{\text{relax}}} \approx \frac{1.2}{1 + N_p} \left( \frac{m}{10 M_{\odot}} \right)^{-1} \left( \frac{M}{10^6 M_{\odot}} \right)^{0.06} . \quad (26) $$

Although Eq. (26) is valid for an isothermal density profile and employed a specific $M-\sigma$ relation [98], we verified that for more sophisticated Dehnen profiles [99] and alternative scaling relations the result holds within a factor of 2. The ratio is roughly independent of mass (but does depend on the adopted $M-\sigma$ relation), and most importantly, it is larger than unity only if $N_p \approx 0$. In this case, the depletion time is longer than the relaxation time and we can therefore assume that the EMRI rate is sustainable. However, since in general there are several plunges per EMRI, $t_d/t_{\text{relax}} < 1$ and a steady state situation where the EMRI rate is given by Eq. (23) cannot be sustained. We therefore define a duty cycle

$$ \Gamma = \min \left\{ \frac{t_d}{t_{\text{relax}}}, 1 \right\}, \quad (27) $$

and a sustainable EMRI rate is given by $\Gamma R_0$.

We can now compute a MBH mass growth by combining this rate with the amount of time a given MBH is surrounded by a cusp, and is therefore a potential EMRI host. We define this time as

$$ t_{\text{EMRI}} = \int dz \frac{dt}{dz} p_0(M, z), \quad (28) $$

where $p_0(M, z)$ is given by Eq. (21) and represents the probability that a MBH of a particular mass is hosted in a stellar cusp as a function of redshift. This time is plotted in the lower panel of Fig. 4 and, as expected, is essentially the Hubble time $T_H$ at $M < 10^5 M_{\odot}$ and rapidly drops to 2 Gyr at $M > 10^6 M_{\odot}$. The mass growth is then

$$ \Delta M = m\Gamma R_0 t_{\text{EMRI}}. \quad (29) $$

Even with the corrective factor of Eq. (27), the supply of COs on inspiralling and plunging orbits can overgrow MBHs. If, for example, $t_{\text{EMRI}} = T_H$ and $m = 10 M_{\odot}$ is the mass of the accreted CO, then from Eq. (29) $\Delta M = m\Gamma R_0 T_H \gtrsim M$ for $M \approx 10^4 M_{\odot}$. For a practical computation of the rate we therefore introduce a damping factor to (arbitrarily) cap the maximum allowed mass growth to be $\exp(-1)M$, so that the MBH can at most grow by an e-fold due to CO accretion in its lifetime. We pick $N_p$ and compute $\Gamma$ from Eq. (26) considering a specific $M-\sigma$ relation, and then calculate, for each MBH mass, $\Delta M$ from Eq. (29). Using this, the damping factor is defined as

<figure><img src="image_6.png" /><figcaption>FIG. 4. Top panel: The adjusted EMRI rate computed according to Eq. (31). The three (central) thick lines assume $N_p = 10$ and correspond to the pessimistic (KormendyHo13, short-dashed orange), fiducial (Gultekin09, solid turquoise) and optimistic (GrahamScott13, long-dashed violet) $M-\sigma$ relations. The two thin turquoise lines show the rates for the fiducial model, but assuming $N_p = 0$ (lower curve) and $N_p = 100$ (upper curve). Lower panel: The average time $t_{\text{EMRI}}$ that a MBH of a given mass is surrounded by a stellar cusp, and is therefore a potential EMRI source, as implicitly defined by Eq. (28). The curves are for the same three different $M-\sigma$ relations in the top panel.</figcaption></figure>

$$ \kappa = \min \left\{ \exp(-1) \frac{M}{\Delta M}, 1 \right\}. \quad (30) $$

Incorporating this, the effective EMRI rate is given by:

$$ R = \kappa \Gamma R_0. \quad (31) $$

Examples of the impact of the factors $\Gamma$ and $\kappa$ on the EMRI rate $R$ are shown in the upper panel of Fig. 4. From Eq. (25), it is clear that the value of the $\Gamma$ factor, and hence $R$, depends critically on $N_p$. We therefore explore three different models featuring $N_p = 0, 10, \text{ and } 100$. Since $\Gamma \approx 1$ for $N_p = 0, N_p$ has the obvious effect of renormalizing the EMRI rate. Equation (26) implies that $R_0$ given by Eq. (23) is close to the supply CO rate allowed by relaxation; if a large fraction of those COs result in direct plunges rather than EMRIs, the EMRI rate must drop accordingly. Different $M-\sigma$ relations alter the slope of the rate as a function of mass because of its influence on the relaxation time and so $\Gamma$. The $\kappa$ factor affects the rates mostly for masses below $\sim 10^5 M_{\odot}$, where overgrowth by CO accretion is easy. This makes our EMRI estimates conservative since it implies a rate suppression. We will see
later that most of LISA detections come from systems with $M > 10^5 M_\odot$, and therefore our results are not severely impacted by the introduction of this damping factor.

## 1. Eccentricity and inclination

To estimate the distribution of EMRI eccentricities at the last stable orbit (LSO), we evolved a large sample of COs from their typical capture radius (which is of the order of 0.01 pc for a $10^6 M_\odot$ central MBH [43]) to final plunge. We found a rather flat eccentricity distribution at plunge in the range $0 < e_p < 0.2$, with a small tail of outliers with higher $e_p$. In the following, we therefore assume a flat distribution in the range $0 < e_p < 0.2$ for simplicity.

Finally, the rates in Eq. (23) have been computed in the nonspinning approximation. As a consequence of frame-dragging effects, the location of the LSO of a test particle orbiting a spinning MBH depends on the inclination of its orbital angular momentum with respect to the MBH spin, $\theta_{inc}$, and on whether the orbit is prograde ($0 \le \theta_{inc} \le \pi/2$) or retrograde ($\pi/2 \le \theta_{inc} \le \pi$). Therefore, the rate $R$ has to be adjusted using a spin-dependent and inclination correcting factor $W(a, \theta_{inc})$ [90]. The correction factor is based on the ratio of the semimajor axis for the spinning case $a_{LSO}^{\text{Kerr}}$ with respect to the nonspinning case $a_{LSO}^{\text{Schw}}$ averaged over the eccentricity, that is [90]:

$$W(a, \theta_{\text{inc}}) = \left\langle \frac{a_{\text{LSO}}^{\text{Kerr}}(e)}{a_{\text{LSO}}^{\text{Schw}}(e)} \right\rangle_e . \quad (32)$$

In practice we use an averaged correcting factor $W(a)$ which is the result of averaging over the orbital inclination

$$W(a) = \langle W(a, \theta_{\text{inc}}) \rangle_{\theta_{\text{inc}}}. \quad (33)$$

In terms of this function, the event rates for EMRIs in the spinning case are related to the nonspinning approximation by

$$R^{\text{Kerr}}(a) = R^{\text{Schw}}[W(a)]^{-0.83}, \quad (34)$$

assuming an old, segregated cusp of COs of mass $10 M_\odot$ around the MBH [83,100].

## D. Putting the pieces together

In summary, the EMRI rate depends on a number of ingredients, as we described above:

(i) The MBH mass function, for which we assume two models: Barausse12 and Gair10.

(ii) The MBH spin distribution, for which we explore three cases: the near-maximally spinning distribution (a98); a flat spin distribution (aflat), and non-spinning MBHs (a0).

(iii) The $M-\sigma$ relation, defining the properties of the stellar distribution surrounding the MBH, the cusp

regrowth time following MBH binary erosion, and the EMRI duty cycle. We consider three relations: Gultekin09, KormendyHo13 and GrahamScott13. We also consider an extra model assuming the Gultekin09 relation but with no cusp erosion.

(iv) The ratio of plunges to EMRIs, assumed to be $N_p = 0, 10$ and $100$.

(v) The characteristic CO mass, for which we consider both $m = 10 M_\odot$ and $m = 30 M_\odot$.

Our default model is based on the self-consistent semi-analytic code for MBH formation and cosmic evolution developed in Barausse [67]. The MBH mass function is therefore the Barausse12, and MBHs are consistently maximally spinning (a98). We use the $M-\sigma$ relation Gultekin09 to compute the cusp regrowth time following MBHB mergers and the EMRI duty cycle. We assume a moderately large number of plunges per EMRI $N_p = 10$ and a characteristic CO mass $m = 10 M_\odot$. Starting from this default set-up, we explore the effect of each single ingredient listed above by varying them individually keeping all the other fixed. We further explore the most optimistic and pessimistic models allowed by all the combinations of the ingredients listed above.

In total, we consider 12 models that we label “Mx” with $x = 1, \dots, 12$. The default setup described above is indicated as M1, and the key to read the models and their main properties are listed in Table I.

For each model we construct the population of EMRIs by Monte Carlo sampling from the distribution $d^3N/(dMdzda) \times p_0(M, z)R(M, a)$. This gives a catalog of EMRIs including the two masses $(M, m)$, redshift of the event $z$, and MBH spin $a$. To define each individual event and construct EMRI waveforms we need to specify 10 more parameters:

(i) Phase, sky position and orientation angles: we assume that the sky position and spin orientation vectors are distributed isotropically on the sphere. The three phases at plunge corresponding to orbital phase, phase of precession of the periapsis and phase of precession of the orbital plane are uniformly distributed between 0 and $2\pi$.

(ii) Inclination and eccentricity are distributed as described in Sec. III C 1.

(iii) Plunge times are taken to be uniform in $[0, 2]$ yr. We ignore events that plunge after the end of the mission duration, although they might be detectable if they are close enough.

Table I illustrates the potential range in the intrinsic EMRI rate. The last column lists the number of EMRIs occurring in the Universe in 1 year (observed at Earth) up to $z = 4.5$ (for model M4 we also report the rate up to $z = 6.5$ in parentheses).

Numbers span more than three order of magnitudes, ranging from about 10 to 20000. The variation is mostly due to the unknown number of plunges and to the poorly
<table><caption>TABLE I. List of EMRI models considered in this work. Column 1 defines the label of each model. For each model we specify the MBH mass function (column 2), the MBH spin model (column 3), whether we consider the effect of cusp erosion following MBH binary mergers (column 4), the *M*–*σ* relation (column 5), the ratio of plunges to EMRIs (column 6), the mass of the COs (column 7); the total number of EMRIs occurring in a year up to *z* = 4.5 (column 8; for model M4 we also show the total rate per year up to *z* = 6.5); the detected EMRI rate per year, with AKK (column 9) and AKS (column 10) waveforms. The AKK and AKS waveforms are introduced in Sec. IV, and bracket waveform modeling uncertainties.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Mass function</th><th rowspan="2">MBH spin</th><th rowspan="2">Cusp erosion</th><th rowspan="2">M-σ relation</th><th rowspan="2">N<sub>p</sub></th><th rowspan="2">CO mass [M<sub>☉</sub>]</th><th colspan="3">EMRI rate [yr<sup>-1</sup>]</th></tr><tr><th>Total</th><th>Detected (AKK)</th><th>Detected (AKS)</th></tr></thead><tbody><tr><td>M1</td><td>Barausse12</td><td>a98</td><td>yes</td><td>Gultekin09</td><td>10</td><td>10</td><td>1600</td><td>294</td><td>189</td></tr><tr><td>M2</td><td>Barausse12</td><td>a98</td><td>yes</td><td>KormendyHo13</td><td>10</td><td>10</td><td>1400</td><td>220</td><td>146</td></tr><tr><td>M3</td><td>Barausse12</td><td>a98</td><td>yes</td><td>GrahamScott13</td><td>10</td><td>10</td><td>2770</td><td>809</td><td>440</td></tr><tr><td>M4</td><td>Barausse12</td><td>a98</td><td>yes</td><td>Gultekin09</td><td>10</td><td>30</td><td>520 (620)</td><td>260</td><td>221</td></tr><tr><td>M5</td><td>Gair10</td><td>a98</td><td>no</td><td>Gultekin09</td><td>10</td><td>10</td><td>140</td><td>47</td><td>15</td></tr><tr><td>M6</td><td>Barausse12</td><td>a98</td><td>no</td><td>Gultekin09</td><td>10</td><td>10</td><td>2080</td><td>479</td><td>261</td></tr><tr><td>M7</td><td>Barausse12</td><td>a98</td><td>yes</td><td>Gultekin09</td><td>0</td><td>10</td><td>15800</td><td>2712</td><td>1765</td></tr><tr><td>M8</td><td>Barausse12</td><td>a98</td><td>yes</td><td>Gultekin09</td><td>100</td><td>10</td><td>180</td><td>35</td><td>24</td></tr><tr><td>M9</td><td>Barausse12</td><td>aflat</td><td>yes</td><td>Gultekin09</td><td>10</td><td>10</td><td>1530</td><td>217</td><td>177</td></tr><tr><td>M10</td><td>Barausse12</td><td>a0</td><td>yes</td><td>Gultekin09</td><td>10</td><td>10</td><td>1520</td><td>188</td><td>188</td></tr><tr><td>M11</td><td>Gair10</td><td>a0</td><td>no</td><td>Gultekin09</td><td>100</td><td>10</td><td>13</td><td>1</td><td>1</td></tr><tr><td>M12</td><td>Barausse12</td><td>a98</td><td>no</td><td>Gultekin09</td><td>0</td><td>10</td><td>20000</td><td>4219</td><td>2279</td></tr></tbody></table>

constrained MBH mass function at $M < 10^6 M_{\odot}$. Cusp erosion has a relative minor effect on the rates (a factor of 2).<fnref n="2" /> Even smaller is the effect of spin, affecting EMRI rates at the 10% level; there are more EMRIs when spins are higher as the LSO is smaller (and so it is more difficult to directly plunge [90]), but this only affects a small portion of orbits. However, we will see that spins will play a more important role in the detectability of these events by LISA. Changing the *M*–*σ* relation, which sets the relation between the MBH and its surrounding population of COs, can introduce a variation of about a factor of 2. More significant are the mass of the COs and the number of plunges, as both of these directly impact the mass accreted by the MBH and so the necessary duty factor to preserve the population of MBHs. An increase in either *m* or *N*<sub>p</sub> by a factor of *X* reduces the EMRI rate by a similar factor. Since we are more uncertain of the number of plunges, this has a greater potential impact on the expected rate, here changing it by almost two orders of magnitude. A drop of about one order of magnitude is achieved by switching to the pessimistic MBH mass distribution, as the reduction in the number of MBHs naturally decreases the number of EMRIs.

For each of the 12 models outlined above we generate 10 Monte-Carlo realizations of the expected population of EMRIs plunging in 1 year. We therefore construct a library of 120 catalogs that includes all EMRI events occurring in the Universe in 10 years for the 12 models.

## IV. WAVEFORMS, SIGNAL ANALYSIS AND PARAMETER ESTIMATION

Having generated astrophysical populations of EMRI systems, we need to determine which of the systems will be observed by LISA. To do this, we need a model of the GW emission from an EMRI system. Accurate gravitational waveforms from EMRIs can be computed using BH perturbation theory, exploiting the large difference in masses of the two objects to regard the smaller as a perturbation of the spacetime of the larger and construct an expansion in the mass ratio (see [101] for a review). Perturbative calculations have not yet been completed to the order necessary to accurately track the phase of an EMRI over an entire inspiral, and these calculations are extremely computationally expensive. Two approximate EMRI models have therefore been developed, which capture the main features of EMRI waveforms at much lower computational cost and can therefore be used to explore questions connected to the detection and scientific exploitation of EMRI observations. Of the two models, the numerical kludge [102,103] is the more accurate and is based on modeling the trajectory of the smaller object as a geodesic of the Kerr background, with inspiral imposed on the system. With further enhancements, the numerical-kludge model may be accurate enough for use in LISA data analysis. However, it is still relatively computationally expensive. The analytic kludge (AK) model [51] is computationally cheaper, at the cost of less faithfulness to real EMRI signals. The AK model approximates gravitational wave emission by that from a Keplerian orbit [104], with precession of the orbital perihelion, precession of the orbital plane, and inspiral of the orbit added using post-Newtonian prescriptions. The AK model provides only an

<footnote n="2">²This could be up to a factor of 4 if kick velocities of few hundred km s⁻¹ are considered in the computation of the cusp regrowth time scale [cf. Eq. (9)].</footnote>
approximation to the true strong-field dynamics, as the orbital frequencies do not exactly match [105], and features such as the final plunge or transient resonances [64,106] are neglected and cannot be readily incorporated. However, the model is cheap to generate and it should include the most important qualitative features of real EMRI signals. The simplicity of the model allows it to be generated in the large numbers required to examine EMRI science questions such as those being explored in this paper, and so we use it here. The AK model has been widely used for similar applications in the literature, in particular it was the EMRI model used in the context of the Mock LISA Data Challenges (MLDCs) [107–110].

The AK model is known to be imperfect, and so in order to quantify inaccuracies we consider two different variants. In the classic work by Barack and Cutler [51], the AK model was cut off when the orbital frequency reached the value corresponding to the Schwarzschild LSO. We denote this form of the AK model by “AKS”, where the “S” stands for “Schwarzschild.” Prograde inspirals into spinning MBH can get much closer before plunge, generating many cycles of higher frequency and amplitude. Thus, omitting those cycles from the model is likely to significantly underestimate the possible signal-to-noise ratio (SNR). An alternative is to continue the inspiral until the frequency reaches the Kerr LSO. We denote this form of the AK model by “AKK”, where the “K” stands for “Kerr.” The post-Newtonian evolution equations used to construct the AK model are increasingly inaccurate as the orbital separation decreases, and so the additional portion of inspiral included in the AKK model is unlikely to be accurately represented, and most likely will lead to an overestimate of the SNR. We will present results for both the AKK model and the AKS model in order to quantify the uncertainty that comes from the modeling assumptions. SNRs can also be computed using results from BH perturbation theory, in particular solutions to the Teukolsky equation, which provides the first-order radiative part of the perturbative evolution. Teukolsky results for circular, equatorial inspirals into spinning BHs were presented in Finn and Thorne [111], and we can use those results to assess the accuracy of the AKS and AKK prescriptions.

Finn and Thorne [111] tabulate their results in terms of corrections relative to a Newtonian inspiral. By setting those corrections equal to 1 we can obtain SNRs for Newtonian inspirals, which we can terminate at the Schwarzschild LSO or at the Kerr LSO. This provides an approximation to the AK model, which is built on Newtonian inspirals, albeit with precession added and inspiral augmented by higher order corrections. Figure 5 shows the sky-averaged horizon distance for a prograde, circular, equatorial inspiral into a black hole with spin $a = 0.99$, computed either using the Teukolsky fluxes, or using Newtonian inspirals truncated at the two different

<figure><img src="image_7.png" /><figcaption>FIG. 5. The redshift at which the sky-averaged SNR of a prograde, circular, equatorial EMRI into a MBH with spin $a = 0.99$ reaches the threshold $q = 20$. The horizon is shown as a function of intrinsic MBH mass and for the two different choices of the compact object mass used in these studies, $m = 10 M_{\odot}$ and $m = 30 M_{\odot}$. The horizon is computed using accurate Teukolsky fluxes and using a Newtonian inspiral truncated either at the Schwarzschild LSO, labeled “AKS,” or at the Kerr LSO, labeled “AKK.” Individual sources may be detected to even larger distances if their orientation is near optimal.</figcaption></figure>

LSOs. We see that, as expected, the approximate AKS and AKK horizons bracket the accurate Teukolsky horizon. The AKS horizon suggests increased sensitivity to lower mass black holes, while the AKK horizon has peak sensitivity at the same MBH mass as the Teukolsky horizon. Although these are just approximations to the true AKS and AKK horizons, we expect the true horizons to have the same shape with the AKS horizon extending to slightly higher redshift than the Newtonian calculations indicate and the AKK horizon to slightly lower redshift, still bracketing the true horizon.

Given a waveform model, we represent the sensitivity of LISA to a given EMRI by a simple SNR threshold. If the EMRI has SNR above the specified threshold, the system will be detected, otherwise it will not. Early work on EMRIs assumed that an SNR of 30 would be required for detection, to allow for the complexities of LISA data analysis [62]. However, in the Mock LISA Data Challenges EMRI signals with SNRs as low as ~15 were successfully identified, albeit under idealized conditions [110]. Therefore, we use a more modest SNR threshold of 20. The SNR is calculated as

$$q = \langle h|h \rangle^{1/2} \qquad (35)$$

using the noise-weighted inner product [112]

$$\langle g|h \rangle = 2 \int_{0}^{\infty} \frac{\tilde{g}(f) \tilde{h}^{*}(f) + \tilde{g}^{*}(f) \tilde{h}(f)}{S_{n}(f)} \mathrm{d}f, \quad (36)$$

where the EMRI signal is denoted by $h(t; \mathbf{\Theta})$, $\mathbf{\Theta}$ represents the parameters of the signal, a tilde indicates the Fourier
transform of the signal, and $S_n(f)$ is the noise power spectral density.

In the limit of suitably high SNR [113], the likelihood for the parameters can be approximated as a Gaussian [114]

$$ \mathcal{L}(\Theta) \propto \exp \left( -\frac{1}{2} \sum_{i,j} \left\langle \frac{dh}{d\Theta^i} \middle| \frac{dh}{d\Theta^j} \right\rangle \Delta\Theta^i \Delta\Theta^j \right), \quad (37) $$

where $\Delta\Theta^i$ represents the displacement from the peak of the distribution (which coincides with the true value in this approximation) for the $i$th parameter. The Fisher matrix has elements

$$ \Gamma_{ij} = \left\langle \frac{dh}{d\Theta^i} \middle| \frac{dh}{d\Theta^j} \right\rangle, \quad (38) $$

and the covariance matrix (the inverse of the Fisher matrix) gives the Cramér-Rao bound on the true width of the distribution [113]. The variance (uncertainty squared) for the $i$th parameter can be approximated by $\sigma_i^2 = (\Gamma^{-1})_{ii}$.

## V. RESULTS

With a number of astrophysically motivated EMRI populations and models of the EMRI waveforms in hand, we are all set to investigate the performance of the LISA detector. The AKS and AKK waveforms introduced in Sec. IV are likely to respectively under- and overestimate

typical EMRI SNRs. In the absence of a more accurate, computationally inexpensive waveform model, we present results for both, with the understanding that they likely bracket the true performance.

### A. Detection rates

To convert from the intrinsic number of EMRIs summarized in Table I to the number of LISA detections, we must compute the SNR of the GW signal and compare to a detection threshold, which we take to be an SNR of $\rho = 20$, as mentioned earlier. The SNR calculation depends on the waveform model; we expect the AKK waveforms to produce larger SNRs due to the extrapolation to the Kerr LSO, resulting in more detectable EMRIs and up to higher redshifts. The SNR distribution of the events above detection threshold is shown in Fig. 6. The distribution approximately follows the characteristic $dN/d\log\rho \propto \rho^{-3}$ behavior of sources uniformly distributed in (Euclidean) volume [115], with small deviations due to cosmological evolution. As expected AKK waveforms lead to a larger number of detections, which is reflected in the higher normalization of the distribution. For models predicting several hundred sources per year (cf. Table I), we predict few events in the tail of the distribution, extending to SNR of a few hundred.

Figure 7 shows the number of detectable signals by using the AKS and AKK waveform models (the rates are also reported in the last two columns of Table I), and compares

<figure><img src="image_8.png" /><figcaption>FIG. 6. SNR distribution for detectable events with AKS and AKK waveforms for all considered models.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Top panel: Event rates (detected and intrinsic) as function of the astrophysical model. Lower panel: Fraction of detection with $\rho > 20$ with respect to the total number of EMRIs featuring a central MBH with $10^4 M_{\odot} < M < 10^7 M_{\odot}$, considering all events at $z < 4.5$ ($z < 6.5$ for M4).</figcaption></figure>

these to the total intrinsic rates reported in Table I. Based on the Teukolsky horizons shown in Fig. 5, and to save computational time, we considered EMRI populations up to $z = 4.5$ for all models with CO mass of $10 M_{\odot}$, and up to $z = 6.5$ when the CO mass is $30 M_{\odot}$ (M4). We will see below that these maximum redshifts are not sufficient to capture all systems detectable using AKK waveforms. We consider this acceptable since AKK waveform generally overestimate EMRI SNR, and the number of missing events

amount to at most a few percent, and thus do not significantly impact our results. As expected, the rates calculated with the AKK model are generally larger because they produce larger SNRs for spinning MBHs. Models M10 and M11 predict the same detectable rates with AKS and AKK waveforms, since they assume that the MBH spins are zero, in which case AKS and AKK waveforms coincide. When using the AKS model the fraction of detectable events is about 10%, independent on the exact features of the model, except for M4 where it increases to around 35%. For the AKK waveform, different spin distributions result in different detection fractions, but these still fall between 10% and 20% in most cases. The expected detection rate is therefore roughly proportional to the intrinsic EMRI rate and ranges between $1 \text{ yr}^{-1}$ and $2000 \text{ yr}^{-1}$ due to severe uncertainties in EMRI astrophysics and dynamics, as discussed in Sec. III.

The fraction of detected events provides clear evidence that EMRI distributions are largely self-similar across the different models, which is confirmed by the (source-frame) mass and redshift distributions of the detected events shown in Fig. 8 and Fig. 9. The sharp $z = 4.5$ cut-off for the AKK case is due to the maximum redshift of the generated population and not to an intrinsic limitation in the detectability of high redshift sources; the small fraction of the number of missing sources should not significantly impact our results. The most common MBH mass is typically between $10^5 M_{\odot}$ and $10^6 M_{\odot}$ in all models. The results

<figure><img src="image_10.png" /><figcaption>FIG. 8. MBH mass distribution for detectable events with AKS and AKK waveforms.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Redshift distribution for detectable events with AKS and AKK waveforms. A maximum source redshift of 4.5 is assumed for all models except M4, where the maximum redshift is 6.5.</figcaption></figure>

based on the AKK waveforms show the detection of more EMRIs into MBHs of larger mass (up to $10^7 M_\odot$), when MBHs are spinning. This is because for such high mass MBHs a prograde inspiral generates a significant number of waveform cycles between the Schwarzschild LSO frequency and the final plunge, and these cycles are at frequencies in the most sensitive range for the LISA detector. Thus, the AKS waveforms omit a significant fraction of the SNR for such systems and underestimates their detectability. This extra contribution to the SNR also allows sources to be seen to further redshift, as illustrated in Fig. 9.

Taken together, Figs. 8 and 9 show that EMRI observations will cover MBHs of $3 \times 10^4 M_\odot < M < 3 \times 10^6 M_\odot$ over a redshift range that is broadly peaked at $0.5 < z < 2$, thus probing a region of the MBH mass–redshift plane that is complementary to both electromagnetic probes of galactic nuclei and LISA observations of MBH binaries. Conventional electromagnetic observations at these low masses out to $z \approx 2$ are extremely challenging, whereas the bulk of LISA MBH binary observations are expected to be at $z > 5$, with only few events expected at $z < 2$ (cf. [116]). EMRIs are a unique opportunity to obtain a large sample of confirmed MBHs at relatively low redshift. Figure 8 further highlights that the number of detected EMRIs is sensitive to the minimum mass scale of nuclear MBHs (Alexander and Bar-Or [93] recently proposed a

universal lower limit of about $2 \times 10^5 M_\odot$), but in the majority of the investigated models, we predict a few detections at $M > 10^6 M_\odot$, which is a relatively safe mass range as it has already been explored by MBH measurements in the local Universe (see, e.g., [117]).

Examples of LISA's completeness as an EMRI survey are given in Fig. 10, where we plot the fraction of detected sources in the (source-frame) mass–redshift plane for selected models. In the default M1 case, LISA will provide an essentially complete survey in the $10^5 M_\odot–10^6 M_\odot$ mass range, out to $z \approx 1$, and it is still 50% complete at $z \approx 3$ when AKK waveforms are considered. If inspiralling COs are massive (M4), the survey is complete out to $z \approx 2$ and still 50% complete out beyond $z \approx 4$ for AKK waveforms.

## B. Parameter estimation

Typical EMRIs spend $O(10^5)$ orbits in the LISA frequency band, and key parameters of the system are encoded in the fine details of the waveform phasing modulation (see, e.g., [51]). The redshifted MBH mass $M_z$ sets the characteristic observed frequencies for the EMRI. The rate of inspiral is controlled by the mass ratio, and so gives constraints on the redshifted CO mass $m_z$. The MBH spin $a$ also influences the orbital frequencies, and becomes more important as the inspiral gets closer to the MBH; the spin sets the LSO and the transition to
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10. Fraction of the intrinsic EMRI population detectable by LISA as a function of source-frame total mass and redshift, for models M1 and M4 and with AKS and AKK waveforms.</figcaption></figure>

plunging.<fnref n="3" /> The orbital eccentricity *e* also affects the orbital frequencies; GW emissions tend to circularize the orbit, so eccentricity is more noticeable earlier in the inspiral.

The large number of cycles completed during the inspiral allow us to obtain exquisite constraints on all the intrinsic parameters, as shown in Fig. 11.<fnref n="4" /> Even in the conservative AKS case, the median relative error on both redshifted masses is in the range 10⁻<fnref n="4" />–10⁻⁵ for essentially all models; the spin of the central MBH and the eccentricity at plunge are measured to an absolute precision of about 10⁻<fnref n="4" /> and 10⁻⁵ respectively. Parameter-estimation precisions for the intrinsic parameters are generally better when calculated using AKK waveforms than the AKS waveforms (except for the nonspinning models M10 and M11, for which the AKS and AKK waveforms, and their parameter-estimation errors, coincide). This is because of

the additional information coming from the late inspiral near the Kerr LSO. The difference is most pronounced for the CO mass and MBH spin (for which the measurement improves by a factor of ≈30 on average). The difference is less striking for the MBH mass and (especially) the eccentricity: estimation for the latter improves on average by less than a factor of 10 because it is mostly constrained by the early inspiral.

Extrinsic parameters such as sky location and distance are primarily determined through the signal amplitude and its modulation as LISA orbits the Sun. These parameters are not strongly dependent upon the GW phase, and hence the large number of cycles completed by an EMRI does not translate to high-precision measurements here. Figure 12 shows that the precision of extrinsic parameter measurement is essentially insensitive to the waveform model. The SNR of an individual source may be higher using the AKK waveform, but the overall distribution of SNRs remains largely the same as more quiet signals become detectable [115], and the typical precision in parameter determination is unaffected. On average, the luminosity distance is measured to 5%–10% precision. The luminosity distance is required to convert the observed redshifted masses back to their true source values. Distance uncertainty will therefore be the dominant source of uncertainty in mass measurements.

<footnote n="3">³The end of the waveform, when the CO plunges into the MBH, also encodes some information which is not captured by these Fisher-matrix estimates [118]. However, the instantaneous SNR in an EMRI is sufficiently low that the plunge is not well resolved, and therefore the inclusion of the plunge should make little difference to parameter estimation.</footnote>

<footnote n="4">⁴All uncertainties are 1σ values, except for the sky-localization error ΔΩ, which is the area of an error ellipse for which there is a probability exp(−1) of the source being outside of it.</footnote>
<figure><img src="image_14.png" /><figcaption>FIG. 11. Distribution over observed EMRIs of the expected statistical errors (1σ uncertainties as computed using the Fisher matrix) in the measurement of intrinsic parameters: central MBH redshifted mass (top left), spin (top right), CO mass (bottom left) and eccentricity at plunge (bottom right). The dashed lines mark the first, second and third quartile of the distributions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 12. Distribution of the statistical errors in the measurement of EMRI extrinsic parameters: luminosity distance (left panel; 1σ uncertainty as computed using the Fisher matrix) and sky localization [right panel; the area of an ellipse with probability 1 − exp(−1) of containing the source]. The dashed lines mark the first, second and third quartile of the distributions. In the plot for the sky position, a horizontal solid red line marks an error of 10 deg$^2$.</figcaption></figure>
The sky localization is usually better than 10 deg². This is the typical field of view of future large optical and radio facilities such as the Large Synoptic Survey Telescope [119] or the Square Kilometre Array [120]. EMRIs localized to this accuracy can therefore be covered with a single pointing to check for the possible presence of electromagnetic counterparts, which could be associated with the interaction between the CO and an MBH accretion disk [56,58–60]. Electromagnetic counterparts would be easiest to observe from close by sources, which would also be the loudest, and so the best localized (usually to better than 1 deg²). Identifying a source galaxy from an electromagnetic counterpart would allow for an independent redshift measurement, which would improve the precision of the (source-frame) mass measurements.

Finally, the precise measurements provided by EMRI observations allow us to maps the spacetime of the MBH and check its Kerr nature. The multipolar structure of the Kerr metric is completely determined by its mass and spin. According the no-hair theorem, the quadrupole moment is given by $Q_K = -a^2 M^3$ [121] (see e.g. [54] for a review of tests of the no-hair theorem with LISA). Since EMRIs are expected to probe the multipolar structure of the central MBH spacetime to high accuracy [46–50], they will be able to confirm if the quadrupole moment obeys the expected Kerr relation [122]. In Fig. 13 we show the precision with which possible deviations of $Q$ away from the Kerr quadrupole can be constrained. We plot the error on the dimensionless quantity $Q = (Q - Q_K)/M^3$ (which is independent of the redshifting of masses). We do not consider any particular modified theory of gravity: the parameter $Q$ is just a phenomenological parametrization of hypothetical deviations from the general-relativistic quadrupole moment, and we are interested in determining

<figure><img src="image_16.png" /><figcaption>FIG. 13. Distribution of the statistical error (1σ uncertainties as computed using the Fisher matrix) in measurement of the deviation of the MBH's quadrupole moment away from the Kerr value. The dashed lines mark the first, second and third quartile of each distribution.</figcaption></figure>

what level of deviation would be measurable. As expected, $Q$ is better constrained by using AKK waveforms, since the effect of a modified quadrupole become important only at small distances from the MBH, i.e. in the late inspiral and plunge.

Overall, for all the parameters that we considered, the distributions of the errors are broadly consistent between the different population models. The populations control the number of events, and so are important for considering how much we could learn about the population of MBHs and their host environments, but do not have a significant impact on our ability to extract the parameters for individual EMRIs.

## VI. CONCLUSIONS

In this paper we have performed a comprehensive analysis of the performance of the recently proposed LISA mission with regards to the detection and parameter estimation of EMRIs. For the first time we have attempted to thoroughly investigate the astrophysical uncertainties that affect the calculations of the expected intrinsic EMRI rate. In more detail, we have constructed competing astrophysical models for the EMRI rate as a function of cosmic time, accounting for: the uncertainty on the expected MBH spin magnitude; the disruption of stellar cusps due to mergers; the MBH growth due to EMRIs and plunges of stellar-mass CO's; and possible viable competing choices for the MBH mass function, the CO mass, and the correlation between MBH masses and stellar velocity dispersions. Although simple, our models capture the diversity of plausible astrophysical uncertainties. Overall, we find that these astrophysical assumptions produce a variance of up to three orders of magnitude in the expected intrinsic EMRI rate.

For each astrophysical model, we have computed the number of expected detections with the LISA interferometer, as well as the precision with which the source parameters (both intrinsic and extrinsic) can be recovered. To this purpose, because of computational-time limitations, we have used two time-inexpensive kludge waveform models [51] that we expect should bracket the results that would be obtained with more sophisticated Teukolsky or self-force based templates (cf. Fig. 5). Our main findings are:

(1) Irrespective of the astrophysical model, at least a few EMRIs per year should be detectable by LISA. This number may reach a few thousands per year under the most optimistic astrophysical assumptions.

(2) Except for the most pessimistic astrophysical models, we predict at least a few events per year should be observable with SNR of several hundreds.

(3) The typical (source-frame) mass and redshift range of detected EMRIs will be $M \sim 10^5-10^6 M_\odot$ and $z \lesssim 2-3$, although we may have events with masses an order of magnitude outside of this range or with
larger redshifts (up to $z \sim 4$ and $z \sim 6$ for COs of 10 $M_{\odot}$ and 30 $M_{\odot}$, respectively) in all but the most pessimistic astrophysical models.

(4) Typical fractional statistical errors with which the intrinsic EMRI parameters (redshifted masses, MBH spin and orbital eccentricity) are expected to be recovered are of the order of $10^{-6}-10^{-4}$. Tests of the multipolar structure of the Kerr metric, which only depend upon these mass and spin measurements, can be performed to percent level precision or better. To convert the redshifted masses to the intrinsic source-frame masses requires the luminosity distance, which is typically inferred to 10% precision. Sky localization is usually of the order of a few square degrees. It is crucial to model the gravitational waveforms in the late inspiral near the plunge to accurately extract the intrinsic parameters, but this has little impact on the extrinsic parameters.

These observations could have impact in three distinct areas: astrophysics, cosmology and fundamental physics.

We have seen that LISA will provide precise measurements of the parameters of individual systems, but more information about the astrophysics of these sources will come from studies of populations. It was shown in Gair *et al.* [61] that the observation of just 10 EMRIs with the classic 5 Gm LISA configuration would be sufficient to measure the slope of the MBH mass function in the local Universe to a precision of $\pm 0.3$. This is the level to which it is currently constrained by electromagnetic observations [61]. The precision with which LISA can measure EMRI parameters does not depend strongly on the configuration of the instrument, so this conclusion should carry over to the current analysis. In all the models except the most pessimistic ones, we expect to see many more than 10 EMRIs, so we would expect to be able to do a high precision measurement of the MBH mass function. One caveat is that what we can measure is the convolution of the MBH mass function with the rate of EMRIs per MBH, not the mass function itself. In Gair *et al.* [61] it was assumed that the latter was known, but as we have described here there are many significant uncertainties. It is an open question as to whether these uncertainties can be reduced or at least quantified, or whether LISA observations will be able to decouple them, for instance by using information from the observed MBH mergers. In addition to the MBH mass function, EMRI observations will provide information on the MBH spin distribution, on the properties of the stellar populations in the centers of galaxies and on the relative efficiency of the mechanisms that lead to EMRI formation.

Observations of GW sources provide measurements of the luminosity distance that can be used to measure the expansion history of the Universe [123]. Individual events do not provide redshifts, but such constraints can be determined statistically. In MacLeod and Hogan [124] it was shown that if LISA observed $\sim 20$ EMRI events at a

redshift $z < 0.5$ it would be possible to determine the Hubble constant to better than 1% by using statistical redshifts estimated from galaxy surveys. We find that all but four of our models predict more than 20 EMRI events at $z < 0.5$.<fnref n="5" /> However, in MacLeod and Hogan [124] it was assumed that LISA would determine the luminosity distance and sky location of an EMRI at redshift $z$ to precisions $\Delta(\ln D_L) < 0.07z$ and $\Delta\Omega < 16z^2$, which were appropriate for the classic 5 Gm LISA configuration, but are optimistic for the current configuration [29]. We find that in the models which have 20 EMRIs at $z < 0.5$, there are at least 5 that also meet the assumed error constraint. If we used only the events at $z < 0.5$, and with errors smaller than these bounds, we would therefore expect to determine the Hubble constant to at least $\sim 2\%$. However, the events with larger errors and events at higher redshift will also contribute to the bound, so we are likely to do better than this, and this should be further explored. In addition, our results show that EMRIs could be detected to much higher redshift than once assumed, which will provide constraints on other cosmological parameters.

The final scientific application of EMRI observations is to test fundamental physics. We have already discussed one such application of EMRI observations, the measurement of the quadrupole deviation from the Kerr metric characterized by $Q$. Every EMRI will provide a percent or better constraint on that parameter, which is comparable to the expectations for the classic LISA mission configuration. This is no surprise, as the key requirement for a test of fundamental physics is to track the phase of an EMRI over a full inspiral, which has to be done in order to find the EMRI in the data using matched filtering. Thus, any EMRI that is detected will provide a powerful test of fundamental physics, and all of the tests previously discussed in the literature should be possible (see Gair *et al.* [54] for a review). Our ability to do this science will not be significantly influenced by the particular astrophysical model, although the models that predict larger rates of EMRI events will more likely lead to the detection of a golden EMRI which is particularly close, has high SNR and, therefore, provides particularly strong constraints.

## ACKNOWLEDGMENTS

A. Sesana is supported by the Royal Society. E. Barausse, E. Berti and A. Klein acknowledge support from the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. This work has made use of the Horizon Cluster, hosted by the

<footnote n="5">⁵Models M4, M5, M8 and M11 do not predict more than 20 EMRIs at $z < 0.5$. M5, M8 and (in particular) M11 are the pessimistic models which have few EMRI events overall and are almost certainly conservative. M4 is the model with high mass COs, so although this does not predict many EMRIs at low redshift, the reach is much greater, so there is still strong potential for constraints on other cosmological parameters.</footnote>
Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us. E. Barausse was supported by the APACHE grant (ANR-16-CE31-0001) of the French Agence Nationale de la Recherche. E. Berti was supported by NSF Grant No. PHY-1607130 and by Fundação para a Ciência e a Tecnologia (FCT) Contract No. IF/00797/2014/CP1214/CT0012 under the IF2014 Programme. C. F. Sopuerta acknowledges support from Contracts No. ESP2013-

47637-P and No. ESP2015-67234-P (Spanish Ministry of Economy and Competitivity of Spain, MINECO). C. P. L. Berry is supported by the Science and Technology Facilities Council. P. Amaro-Seoane acknowledges support from the Ramón y Cajal Programme of the Ministry of Economy, Industry and Competitiveness of Spain. This work has been partially supported by the CAS President's International Fellowship Initiative. This work was supported by the Centre National d'Études Spatiales.

[1] J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso et al., *Classical Quantum Gravity* **32**, 074001 (2015).

[2] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. Lett.* **116**, 061102 (2016).

[3] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. Lett.* **116**, 241103 (2016).

[4] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. D* **93**, 122003 (2016).

[5] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. X* **6**, 041015 (2016).

[6] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. Lett.* **116**, 241102 (2016).

[7] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. X* **6**, 041014 (2016).

[8] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Astrophys. J. Lett.* **818**, L22 (2016).

[9] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Phys. Rev. Lett.* **116**, 221101 (2016).

[10] N. Yunes, K. Yagi, and F. Pretorius, *Phys. Rev. D* **94**, 084002 (2016).

[11] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin et al., *Classical Quantum Gravity* **32**, 024001 (2015).

[12] K. Somiya, *Classical Quantum Gravity* **29**, 124007 (2012).

[13] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari et al., *Living Rev. Relativ.* **19** (2016).

[14] D. Lynden-Bell and M. J. Rees, *Mon. Not. R. Astron. Soc.* **152**, 461 (1971).

[15] A. Soltan, *Mon. Not. R. Astron. Soc.* **200**, 115 (1982).

[16] R. Genzel, F. Eisenhauer, and S. Gillessen, *Rev. Mod. Phys.* **82**, 3121 (2010).

[17] J. Kormendy and D. Richstone, *Annu. Rev. Astron. Astrophys.* **33**, 581 (1995).

[18] K. Gültekin, D. O. Richstone, K. Gebhardt, T. R. Lauer, S. Tremaine, M. C. Aller, R. Bender, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, L. C. Ho, J. Kormendy, J. Magorrian, J. Pinkney, and C. Siopis, *Astrophys. J.* **698**, 198 (2009).

[19] J. Kormendy and L. C. Ho, *Annu. Rev. Astron. Astrophys.* **51**, 511 (2013).

[20] A. King, *Astrophys. J. Lett.* **596**, L27 (2003).

[21] Z. Haiman and E. Quataert, in *Supermassive Black Holes in the Distant Universe*, Astrophysics and Space Science Library, edited by A. J. Barger (Springer Netherlands, Dordrecht (The Netherlands), 2004), Vol. 308, p. 147.

[22] M. Volonteri and P. Natarajan, *Mon. Not. R. Astron. Soc.* **400**, 1911 (2009).

[23] R. Schödel, A. Feldmeier, N. Neumayer, L. Meyer, and S. Yelda, *Classical Quantum Gravity* **31**, 244007 (2014).

[24] T. Alexander, *Phys. Rep.* **419**, 65 (2005).

[25] J. G. Hills, *Astrophys. J.* **80**, 809 (1975).

[26] B. W. Murphy, H. N. Cohn, and R. H. Durisen, *Astrophys. J.* **370**, 60 (1991).

[27] M. Freitag and W. Benz, *Astron. Astrophys.* **394**, 345 (2002).

[28] S. Gezari, J. P. Halpern, S. Komossa, D. Grupe, and K. M. Leighly, *Astrophys. J.* **592**, 42 (2003).

[29] *Amaro-Seoane et al.*, arXiv:1702.00786.

[30] M. Armano, H. Audley, G. Auger, J. T. Baird, M. Bassan, P. Binetruy, M. Born, D. Bortoluzzi, N. Brandt et al., *Phys. Rev. Lett.* **116**, 231101 (2016).

[31] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder, F. Ohme, and B. Wardell, *Phys. Rev. D* **93**, 024003 (2016).

[32] C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller, G. Servant, and D. J. Weir, *J. Cosmol. Astropart. Phys.* **04** (2016) 001.

[33] N. Bartolo, C. Caprini, V. Domcke, D. G. Figueroa, J. Garcia-Bellido, M. Chiara Guzzetti, M. Liguori, S. Matarrese, M. Peloso, A. Petiteau, A. Ricciardone, M. Sakellariadou, L. Sorbo, and G. Tasinato, *J. Cosmol. Astropart. Phys.* **12** (2016) 026.

[34] N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein, and A. Petiteau, *J. Cosmol. Astropart. Phys.* **04** (2016) 002.

[35] E. Barausse, N. Yunes, and K. Chamberlain, *Phys. Rev. Lett.* **116**, 241104 (2016).
[36] E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. **117**, 101102 (2016).

[37] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D **73**, 064030 (2006).

[38] A. Sesana, Phys. Rev. Lett. **116**, 231102 (2016).

[39] A. Nishizawa, E. Berti, A. Klein, and A. Sesana, Phys. Rev. D **94**, 064020 (2016).

[40] A. Nishizawa, A. Sesana, E. Berti, and A. Klein, Mon. Not. R. Astron. Soc. **465**, 4375 (2017).

[41] K. Breivik, C. L. Rodriguez, S. L. Larson, V. Kalogera, and F. A. Rasio, Astrophys. J. Lett. **830**, L18 (2016).

[42] S. Vitale, Phys. Rev. Lett. **117**, 051102 (2016).

[43] P. Amaro-Seoane, J. R. Gair, M. Freitag, M. C. Miller, I. Mandel, C. J. Cutler, and S. Babak, Classical Quantum Gravity **24**, R113 (2007).

[44] P. C. Peters, Phys. Rev. **136**, B1224 (1964).

[45] T. Hinderer and É. É. Flanagan, Phys. Rev. D **78**, 064028 (2008).

[46] F. D. Ryan, Phys. Rev. D **52**, 5707 (1995).

[47] F. D. Ryan, Phys. Rev. D **56**, 1845 (1997).

[48] F. D. Ryan, Phys. Rev. D **56**, 7732 (1997).

[49] K. Glampedakis and S. Babak, Classical Quantum Gravity **23**, 4167 (2006).

[50] L. Barack and C. Cutler, Phys. Rev. D **75**, 042003 (2007).

[51] L. Barack and C. Cutler, Phys. Rev. D **69**, 082005 (2004).

[52] E. A. Huerta and J. R. Gair, Phys. Rev. D **79**, 084021 (2009).

[53] K. G. Arun, S. Babak, E. Berti, N. Cornish, C. Cutler, J. Gair, S. A. Hughes, B. R. Iyer, R. N. Lang, I. Mandel, E. K. Porter, B. S. Sathyaprakash, S. Sinha, A. M. Sintes, M. Trias, C. Van Den Broeck, and M. Volonteri, Classical Quantum Gravity **26**, 094027 (2009).

[54] J. R. Gair, M. Vallisneri, S. L. Larson, and J. G. Baker, Living Rev. Relativ. **16** (2013).

[55] E. Barausse, L. Rezzolla, D. Petroff, and M. Ansorg, Phys. Rev. D **75**, 064026 (2007).

[56] E. Barausse and L. Rezzolla, Phys. Rev. D **77**, 104027 (2008).

[57] J. R. Gair, É. É. Flanagan, S. Drasco, T. Hinderer, and S. Babak, Phys. Rev. D **83**, 044037 (2011).

[58] N. Yunes, B. Kocsis, A. Loeb, and Z. Haiman, Phys. Rev. Lett. **107**, 171103 (2011).

[59] E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D **89**, 104059 (2014).

[60] E. Barausse, V. Cardoso, and P. Pani, J. Phys. Conf. Ser. **610**, 012044 (2015).

[61] J. R. Gair, C. Tang, and M. Volonteri, Phys. Rev. D **81**, 104014 (2010).

[62] J. R. Gair, L. Barack, T. Creighton, C. Cutler, S. L. Larson, E. S. Phinney, and M. Vallisneri, Classical Quantum Gravity **21**, S1595 (2004).

[63] J. R. Gair, Classical Quantum Gravity **26**, 094034 (2009).

[64] C. P. L. Berry, R. H. Cole, P. Cañizares, and J. R. Gair, Phys. Rev. D **94**, 124042 (2016).

[65] P. Amaro-Seoane et al., Classical Quantum Gravity **29**, 124016 (2012).

[66] P. Amaro-Seoane et al., arXiv:1305.5720.

[67] E. Barausse, Mon. Not. R. Astron. Soc. **423**, 2533 (2012).

[68] A. Sesana, E. Barausse, M. Dotti, and E. M. Rossi, Astrophys. J. **794**, 104 (2014).

[69] F. Antonini, E. Barausse, and J. Silk, Astrophys. J. Lett. **806**, L8 (2015).

[70] F. Antonini, E. Barausse, and J. Silk, Astrophys. J. **812**, 72 (2015).

[71] P. Madau and M. J. Rees, Astrophys. J. Lett. **551**, L27 (2001).

[72] F. Shankar, D. H. Weinberg, and J. Miralda-Escudé, Astrophys. J. **690**, 20 (2009).

[73] F. Shankar, Classical Quantum Gravity **30**, 244001 (2013).

[74] M. Dotti, M. Colpi, S. Pallini, A. Perego, and M. Volonteri, Astrophys. J. **762**, 68 (2013).

[75] J. N. Bahcall and R. A. Wolf, Astrophys. J. **216**, 883 (1977).

[76] M. Milosavljević and D. Merritt, Astrophys. J. **563**, 34 (2001).

[77] J. Binney and S. Tremaine, *Galactic Dynamics* (Princeton University Press, Princeton, NJ, 1987), p. 747.

[78] A. Sesana, F. Haardt, and P. Madau, Astrophys. J. **651**, 392 (2006).

[79] D. C. Heggie, Mon. Not. R. Astron. Soc. **173**, 729 (1975).

[80] J. G. Hills, Astron. J. **80**, 809 (1975).

[81] G. D. Quinlan, New Astron. **1**, 35 (1996).

[82] A. Sesana and F. M. Khan, Mon. Not. R. Astron. Soc. Lett. **454**, L66 (2015).

[83] P. Amaro-Seoane and M. Preto, Classical Quantum Gravity **28**, 094017 (2011).

[84] F. M. Khan, M. Preto, P. Berczik, I. Berentzen, A. Just, and R. Spurzem, Astrophys. J. **749**, 147 (2012).

[85] A. W. Graham and N. Scott, Astrophys. J. **764**, 151 (2013).

[86] F. Shankar, M. Bernardi, R. K. Sheth, L. Ferrarese, A. W. Graham, G. Savorgnan, V. Allevato, A. Marconi, R. Läsker, and A. Lapi, Mon. Not. R. Astron. Soc. **460**, 3119 (2016).

[87] D. Merritt, Rep. Prog. Phys. **69**, 2513 (2006).

[88] C. Hopman and T. Alexander, Astrophys. J. **645**, 1152 (2006).

[89] D. Merritt, T. Alexander, S. Mikkola, and C. M. Will, Phys. Rev. D **84**, 044024 (2011).

[90] P. Amaro-Seoane, C. F. Sopuerta, and M. D. Freitag, Mon. Not. R. Astron. Soc. **429**, 3155 (2013).

[91] D. Merritt, Astrophys. J. **814**, 57 (2015).

[92] B. Bar-Or and T. Alexander, Astrophys. J. **820**, 129 (2016).

[93] T. Alexander and B. Bar-Or, arXiv:1701.00415.

[94] E. E. Salpeter, Astrophys. J. **121**, 161 (1955).

[95] P. Kroupa, Mon. Not. R. Astron. Soc. **322**, 231 (2001).

[96] S. E. Woosley, A. Heger, and T. A. Weaver, Rev. Mod. Phys. **74**, 1015 (2002).

[97] M. Spera, M. Mapelli, and A. Bressan, Mon. Not. R. Astron. Soc. **451**, 4086 (2015).

[98] S. Tremaine, K. Gebhardt, R. Bender, G. Bower, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, C. Grillmair, L. C. Ho, J. Kormendy, T. R. Lauer, J. Magorrian, J. Pinkney, and D. Richstone, Astrophys. J. **574**, 740 (2002).

[99] W. Dehnen, Mon. Not. R. Astron. Soc. **265**, 250 (1993).

[100] T. Alexander and C. Hopman, Astrophys. J. **697**, 1861 (2009).

[101] E. Poisson, A. Pound, and I. Vega, Living Rev. Relativ. **14**, 7 (2011).

[102] J. R. Gair and K. Glampedakis, Phys. Rev. D **73**, 064037 (2006).
[103] S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes, Phys. Rev. D **75**, 024005 (2007).

[104] P. C. Peters and J. Mathews, Phys. Rev. **131**, 435 (1963).

[105] A. J. K. Chua and J. R. Gair, Classical Quantum Gravity **32**, 232002 (2015).

[106] É. É. Flanagan and T. Hinderer, Phys. Rev. Lett. **109**, 071102 (2012).

[107] K. A. Arnaud, S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, C. Cutler, S. L. Larson, B. S. Sathyaprakash, M. Vallisneri, A. Vecchio, and J.-Y. Vinet, in *Laser Interferometer Space Antenna: 6th International LISA Symposium*, American Institute of Physics Conference Series, Vol. 873, edited by S. M. Merkovitz and J. C. Livas (AIP Publishing, Melville (New York, USA), 2006), p. 625.

[108] K. A. Arnaud, S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, C. Cutler, L. S. Finn, S. L. Larson, T. Littenberg, E. K. Porter, M. Vallisneri, A. Vecchio, and J.-Y. Vinet, Classical Quantum Gravity **24**, S551 (2007).

[109] S. Babak et al., Classical Quantum Gravity **25**, 184026 (2008).

[110] S. Babak et al., Classical Quantum Gravity **27**, 084009 (2010).

[111] L. S. Finn and K. S. Thorne, Phys. Rev. D **62**, 124021 (2000).

[112] L. S. Finn, Phys. Rev. D **46**, 5236 (1992).

[113] M. Vallisneri, Phys. Rev. D **77**, 042001 (2008).

[114] C. Cutler and É. E. Flanagan, Phys. Rev. D **49**, 2658 (1994).

[115] B. F. Schutz, Classical Quantum Gravity **28**, 125023 (2011).

[116] A. Sesana, J. Gair, E. Berti, and M. Volonteri, Phys. Rev. D **83**, 044036 (2011).

[117] N. Scott, A. W. Graham, and J. Schombert, Astrophys. J. **768**, 76 (2013).

[118] I. Mandel, C. P. L. Berry, F. Ohme, S. Fairhurst, and W. M. Farr, Classical Quantum Gravity **31**, 155005 (2014).

[119] P. A. Abell, J. Allison, S. F. Anderson, J. R. Andrew, J. R. P. Angel, L. Armus, D. Arnett, S. J. Asztalos, T. S. Axelrod et al., arXiv:0912.0201.

[120] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. W. Lazio, IEEE Proc. **97**, 1482 (2009).

[121] R. O. Hansen, J. Math. Phys. (N.Y.) **15**, 46 (1974).

[122] K. Yagi and L. C. Stein, Classical Quantum Gravity **33**, 054001 (2016).

[123] B. F. Schutz, Nature (London) **323**, 310 (1986).

[124] C. L. MacLeod and C. J. Hogan, Phys. Rev. D **77**, 043512 (2008).