# Gravitational wave searches for ultralight bosons with LIGO and LISA

Richard Brito, Shrobana Ghosh, Enrico Barausse, Emanuele Berti, Vitor

Cardoso, Irina Dvorkin, Antoine Klein, Paolo Pani

► To cite this version:

Richard Brito, Shrobana Ghosh, Enrico Barausse, Emanuele Berti, Vitor Cardoso, et al. Gravitational wave searches for ultralight bosons with LIGO and LISA. Physical Review D, 2017, 96 (6), pp.064050. 10.1103/PhysRevD.96.064050. [hal-01645860]

HAL Id: hal-01645860

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

Submitted on 4 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
# Gravitational wave searches for ultralight bosons with LIGO and LISA

Richard Brito,¹,* Shrobana Ghosh,² Enrico Barausse,³ Emanuele Berti,²,⁴ Vitor Cardoso,⁴,⁵
Irina Dvorkin,³,⁶ Antoine Klein,³ and Paolo Pani⁷,⁴

¹Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, Potsdam-Golm 14476, Germany

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

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

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

⁵Perimeter Institute for Theoretical Physics, 31 Caroline Street North Waterloo, Ontario N2L 2Y5, Canada

⁶Institut Lagrange de Paris (ILP), Sorbonne Universités, 98 bis bd Arago, 75014 Paris, France

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

(Received 23 June 2017; published 27 September 2017)

Ultralight bosons can induce superradiant instabilities in spinning black holes, tapping their rotational energy to trigger the growth of a bosonic condensate. Possible observational imprints of these boson clouds include (i) direct detection of the nearly monochromatic (resolvable or stochastic) gravitational waves emitted by the condensate, and (ii) statistically significant evidence for the formation of “holes” at large spins in the spin versus mass plane (sometimes also referred to as “Regge plane”) of astrophysical black holes. In this work, we focus on the prospects of LISA and LIGO detecting or constraining scalars with mass in the range $m_s \in [10^{-19}, 10^{-15}]$ eV and $m_s \in [10^{-14}, 10^{-11}]$ eV, respectively. Using astrophysical models of black-hole populations calibrated to observations and black-hole perturbation theory calculations of the gravitational emission, we find that, in optimistic scenarios, LIGO could observe a stochastic background of gravitational radiation in the range $m_s \in [2 \times 10^{-13}, 10^{-12}]$ eV, and up to $10^4$ resolvable events in a 4-year search if $m_s \sim 3 \times 10^{-13}$ eV. LISA could observe a stochastic background for boson masses in the range $m_s \in [5 \times 10^{-19}, 5 \times 10^{-16}]$, and up to $\sim 10^3$ resolvable events in a 4-year search if $m_s \sim 10^{-17}$ eV. LISA could further measure spins for black-hole binaries with component masses in the range [$10^3, 10^7] M_\odot$, which is not probed by traditional spin-measurement techniques. A statistical analysis of the spin distribution of these binaries could either rule out scalar fields in the mass range $\sim [4 \times 10^{-18}, 10^{-14}]$ eV, or measure $m_s$ with ten percent accuracy if light scalars in the mass range $\sim [10^{-17}, 10^{-13}]$ eV exist.

DOI: 10.1103/PhysRevD.96.064050

## I. INTRODUCTION

The first gravitational wave (GW) detections by the Laser Interferometric Gravitational-wave Observatory (LIGO) are a historical landmark. GW150914 [1], GW151226 [2], GW170104 [3] and the LVT151012 trigger [4] provided the strongest evidence to date that stellar-mass black holes (BHs) exist and merge [5–9]. In this work we discuss the exciting possibility that LIGO and space-based detectors like LISA [10,11] could revolutionize our understanding of dark matter and of fundamental interactions in the Universe.

Ultralight bosons—such as dark photons, the QCD axion or the axionlike particles predicted by the string axiverse scenario—could be a significant component of dark matter [12–15]. These fields interact very feebly with standard

model particles, but the equivalence principle imposes some universality in the way that they gravitate. Light bosonic fields around spinning black holes trigger super-radiant instabilities, which can be strong enough to have astrophysical implications [16]. Therefore, GW detectors can either probe the existence of new particles beyond the standard model or—in the absence of detections—impose strong constraints on their masses and couplings [17–22].

Superradiance by rotating BHs was first demonstrated with a thought experiment involving particles [16,23]. Penrose imagined a particle falling into a BH and splitting into two particles. If the splitting occurs in the ergoregion, one of the fragmentation products can be in a negative-energy state as seen by an observer at infinity, and therefore the other fragmentation product can escape to infinity with energy larger than the original particle. The corresponding process involving waves amplifies any bosonic wave whose frequency $\omega$ satisfies $0 < \omega < m\Omega_H$, where $m$ is the

<footnote>*richard.brito@aei.mpg.de</footnote>
azimuthal index of the (spheroidal) harmonics used to separate the angular dependence, and $\Omega_H$ is the horizon angular velocity [16,24,25]. The wave is amplified at the expense of the BH's rotational energy. If the wave is trapped—for example, through a confining mechanism like a mirror placed at some finite distance—the amplification process will repeat, destabilizing the system. This creates a “BH bomb” [26,27]. Massive fields are naturally trapped by their own mass, leading to a superradiant instability of the Kerr geometry. The time scales and evolution of BH superradiant instabilities were extensively studied by several authors for massive spin-0 [28–31], spin-1 [21,32–35] and spin-2 fields [36], using both analytic and numerical methods.

For a bosonic field with mass $m_s$, superradiant instabilities are strongest when the Compton wavelength of the massive boson $\hbar/(m_s c)$ is comparable to the Schwarzschild radius $R = 2GM/c^2$, where $M$ is the BH mass. Under these conditions the bosonic field can bind to the BH, forming a “gravitational atom.” Instabilities can produce holes in the BH mass/spin plane (sometimes also called the BH “Regge plane”): for a given boson mass, spinning BHs should not exist when the dimensionless spin $\chi \equiv a/M$ is above an instability window centered around values of order unity of the dimensionless quantity [16,17]

$$ \frac{2GMm_s}{c\hbar} = 1.5 \frac{M}{10^6 M_{\odot}} \frac{m_s c^2}{10^{-16} \text{eV}}. \quad (1) $$

Typical instability windows for selected values of $m_s$ are shown as shaded areas in Fig. 1, which shows the spin versus mass plane. These instability windows are obtained by requiring that the instability acts on timescales shorter than known astrophysical processes such as accretion, i.e., we require that the superradiant instability time scales for scalar field perturbations with $l = m = 1, 2, 3$ are shorter than a typical accretion time scale, here conservatively assumed to be the Salpeter time scale defined below for a typical efficiency $\eta=0.1$ and Eddington rate $f_{\text{Edd}} = 1$ [cf. Eq. (51)].

In Fig. 1, black data points denote electromagnetic estimates of stellar or massive BH spins obtained using either the $K\alpha$ iron line or the continuum fitting method [37,38]. Roughly speaking, massive BH spin measurements probe the existence of instability windows in the mass range $m_s \sim 10^{-19}-10^{-17}$ eV. For stellar-mass BHs, the relevant mass range is $m_s \sim 10^{-12}-10^{-11}$ eV. Red data points are LIGO 90% confidence levels for the spins of the primary and secondary BHs in the three merger events detected so far (GW150914, GW151226 and GW170104 [3,4]). For LIGO BH binaries accretion should not be important. In such case, our choice for the reference time scale $t_s$ is conservative: more accurate and stringent constraints can be imposed by comparing the instability time scale with the Hubble time or with the age of the BHs. On the other hand, in Fig. 1 we do not include the remnant BHs detected by LIGO because the observation time scale

<figure><img src="image_2.png" /><figcaption>FIG. 1. Exclusion regions in the BH mass-spin plane (Regge plane) for a massive scalar field. For each mass $m_s$, the instability threshold is obtained by setting the superradiant instability time scales for $l = m = 1, 2, 3$ equal to a typical accretion time scale, taken to be $\tau = 50$ Myr (see main text for details). Black data points (with error bars) are spin estimates of stellar and massive BHs obtained through the $K\alpha$ or continuum fitting methods [37,38]. Red data points are GW measurements of the primary and secondary BHs from the three LIGO detections (GW150914, GW151226 and GW170104 [3,4]). Blue, green and brown data points are projected LISA measurements under the assumption that there are no light bosons for three different astrophysical black hole population models (popIII, Q3 and Q3-nod from [39]), as discussed in the text. We assume a LISA observation time $T_{\text{obs}} = 1$ yr, and to avoid cluttering we only show events for which LISA spin measurement errors are relatively small ($\Delta\chi/\chi \le 2/3$). The top horizontal line is a frequency scale corresponding to the BH mass, $f \approx \mu/\pi$ with $\mu \sim 0.2/M$ as a reference value.</figcaption></figure>

of the latter is obviously much shorter than the superradiant instability time scale, and therefore post-merger observations cannot be used to place constraints on the boson mass.

Blue, green and brown data points are projected LISA measurements for three different astrophysical black-hole population models (popIII, Q3, Q3-nod) from [39], assuming one year of observation. The main point of Fig. 1 is to highlight one of the most remarkable results of this work: LISA BH spin measurements cover the intermediate mass range (roughly $m_s \sim 10^{-13}-10^{-16}$ eV, with the lower and upper bounds depending on the astrophysical model, and more specifically on the mass of BH seeds in the early Universe), unaccessible to electromagnetic observations of stellar and massive BHs. In other words, LISA’s capability to measure the mass and spin of binary BH components out to cosmological distances.<fnref n="1" /> implies that LISA can also

<footnote n="1">¹We do not study holes in the Regge plane for LIGO because spin magnitude measurements for the binary components are expected to be poor, even with third-generation detectors [40,41], and they overlap in mass with existing EM spin estimates.</footnote>
probe the existence of light bosonic particles in a large mass range that is not accessible by other BH-spin measurement methods. In Sec. VI below we quantify this expectation with a more detailed Bayesian model-selection analysis, showing in addition that (if light bosons exist) LISA could measure their mass with ~10% accuracy.

We note that electromagnetic measurements of black-hole spins also provide constraints on the scalar field masses that partly overlap with constraints derived in this paper. For example, the spin measurements of stellar mass BHs disfavor the existence of a scalar field with masses between roughly $2 \times 10^{-11}$ eV > $m_s$ > $6 \times 10^{-13}$ eV [19] and $4 \times 10^{-17}$ eV > $m_s$ > $5 \times 10^{-20}$ eV for massive BHs. However, GW spin measurements and constraints rely on fewer astrophysical assumptions (e.g. on the accretion disk and its spectrum) than electromagnetic constraints, and are therefore more robust. On the other hand, while electromagnetic observations of stellar mass BHs partly overlap with the GW constraints from LIGO, the electromagnetic observation of massive BHs probe lower scalar field masses that the ones coming from GW observations, and are thus complementary to the constraints that we estimate in this paper. Fig. 1 shows that electromagnetic and GW observations should be considered jointly to build evidence for or against the existence of a scalar field with a given mass.

An even more exciting prospect is the direct detection of the GWs produced by a BH-boson condensate system [19–21]. Through superradiance, energy and angular momentum are extracted from a rotating BH and the number of bosons grows exponentially, producing a bosonic “cloud” at distance ~ħ²(2GMm_s²)⁻¹ from the BH. This nonaxisymmetric cloud creates a time-varying quadrupole moment, leading to long-lasting, monochromatic GWs with frequency determined by the boson mass. Thus, the existence of light bosons can be tested (or constrained) directly with GW detectors.

To estimate the detectability of these signals we need careful estimates of the signal strength and astrophysical models for stellar-mass and massive BH populations. Here we compute the GW signal produced by superradiant instabilities using GW emission models in BH perturbation theory [42], which are expected to provide an excellent approximation for all situations of physical interest [18,34,35]. On the astrophysical side, we adopt the same BH formation models [43] that were used in previous LISA studies [39,44–47]. As shown below, semicoherent searches with LISA (LIGO) could detect individual signals at luminosity distances as large as ~2 Gpc (~200 Mpc) for a boson of mass $10^{-17}(10^{-13})$ eV (compare this with the farthest estimated distance for LIGO BH binary merger detections so far, the $880_{-390}^{+450}$ Mpc of GW170104 [3]).

The plan of the paper is as follows. In Sec. II we outline our calculation of gravitational radiation from bosonic condensates around rotating BHs. In Sec. III

and Sec. IV we present our astrophysical models of massive and stellar-mass BH formation, respectively. Our predictions for rates of boson-condensate GW events detectable by LISA and LIGO, either as resolvable events or as a stochastic background, are given in Sec. V. In Sec. VI we use a Bayesian model selection framework to quantify how LISA spin measurements in BH binary mergers can either exclude certain boson mass ranges by looking at the presence of holes in the Regge plane, or (if bosons exist in the Universe) be used to estimate boson masses. We conclude by summarizing our main results and identifying some promising avenues for future work.

In the following, we use geometrized units $G = c = 1$.

## II. GRAVITATIONAL WAVES FROM BOSONIC CONDENSATES AROUND BLACK HOLES

In general, the development of instabilities must be followed through nonlinear evolutions. Numerical studies of the development of superradiant instabilities are still in their infancy (see e.g. [34,35,48–51]), mainly because of the long instability growth time for scalar perturbations, which makes simulations computationally prohibitive. If we restrict attention to near-vacuum environments, the scalar cloud around the spinning BH can only grow by tapping the BH’s rotational energy. Standard arguments [52] imply that the BH can lose at most 29% of its mass. For the process at hand, it turns out that the cloud can store at most ~10% of the BH’s mass [34,53], therefore the spacetime is described to a good approximation by the Kerr metric, and perturbative calculations are expected to give good estimates of the emitted radiation [18,42]. These expectations were recently validated by nonlinear numerical evolutions in the spin-1 case [34,35], where the instability growth time scale is faster. Reassuringly, these numerical simulations are consistent with qualitative and quantitative predictions from BH perturbation theory [21,32,33]. In summary, a body of analytic and numerical work justifies the use of calculations in BH perturbation theory to estimate the gravitational radiation emitted by bosonic condensates around Kerr BHs. We now turn to a detailed description of this calculation.

### A. Test scalar field on a Kerr background

Neglecting possible self-interaction terms or couplings to other fields, the action describing a real scalar field minimally coupled to gravity is

$$S = \int d^4x \sqrt{-g} \left( \frac{R}{16\pi} - \frac{1}{2}g^{\mu\nu}\Psi_{,\mu}\Psi_{,\nu} - \frac{\mu^2}{2}\Psi^2 \right). \quad (2)$$

Here we defined a parameter

$$\mu = m_s/\hbar, \quad (3)$$
which has dimensions of an inverse mass (in our geometrized units). The field equations derived from this action are $\nabla_{\mu}\nabla^{\mu}\Psi = \mu^{2}\Psi$ and $G^{\mu\nu} = 8\pi T^{\mu\nu}$, with

$$T^{\mu\nu} = \Psi^{\mu}\Psi^{\nu} - \frac{1}{2}g^{\mu\nu}(\Psi_{,\alpha}\Psi^{\alpha} + \mu^{2}\Psi^{2}). \quad (4)$$

In the test-field approximation, where the scalar field propagates on a fixed Kerr background with mass $M$ and spin $J = aM$, the general solution of the Klein-Gordon equation can be written as

$$\Psi = \Re \left[ \int d\omega e^{-i\omega t + i\omega\phi} {}_0S_{\ell m\omega}(\vartheta)\psi_{\ell m\omega}(r) \right], \quad (5)$$

where a sum over harmonic indices $(\ell, m)$ is implicit, and ${}_sY_{\ell m\omega}(\vartheta, \phi) = {}_sS_{\ell m\omega}(\vartheta)e^{im\phi}$ are the spin-weighted spheroidal harmonics of spin weight $s$, which reduce to the scalar spheroidal harmonics for $s = 0$ [54]. The radial and angular functions satisfy the following coupled system of differential equations:

$$\begin{aligned} \mathcal{D}_{\vartheta}[{}_0S] &+ \left[ a^2(\omega^2 - \mu^2)\cos^2\vartheta - \frac{m^2}{\sin^2\vartheta} + \lambda \right]{}_0S = 0, \\ \mathcal{D}_r[\psi] &+ [\omega^2(r^2 + a^2)^2 - 4aMrm\omega + a^2m^2 \\ & \qquad - \Delta(\mu^2r^2 + a^2\omega^2 + \lambda)]\psi = 0, \end{aligned}$$

where for simplicity we omit the $(\ell, m)$ subscripts, $r_{\pm} = M \pm \sqrt{M^2 - a^2}$ denotes the coordinate location of the inner and outer horizons, $\Delta = (r - r_+)(r - r_-)$, $\mathcal{D}_r = \Delta\partial_r(\Delta\partial_r)$, and $\mathcal{D}_{\vartheta} = (\sin\vartheta)^{-1}\partial_{\vartheta}(\sin\vartheta\partial_{\vartheta})$. For $a = 0$, the angular eigenfunctions ${}_0S_{\ell m}(\vartheta)$ reduce to the usual scalar spherical harmonics with eigenvalues $\lambda = \ell(\ell + 1)$.

Imposing appropriate boundary conditions, a solution to the above coupled system can be obtained using, e.g. a continued-fraction method [30,31]. Because of dissipation, this boundary value problem is non-Hermitian. The solutions are generically described by an infinite, discrete set of complex eigenfrequencies [55]

$$\omega_{\ell mn} \equiv \omega = \omega_R + i\omega_I, \quad (6)$$

where $n$ is the overtone number and $\{\omega_R, \omega_I\} \in \mathbb{R}$. In particular, this system admits quasi-bound state solutions which become unstable—i.e., from Eq. (5), have $\omega_I > 0$— for modes satisfying the superradiant condition $\omega_R < m\Omega_H$, with $\Omega_H = a/(2Mr_+)$ [28,31]. For these solutions the eigenfunctions are exponentially suppressed at spatial infinity:

$$\psi(r) \propto \frac{r^{\nu} e^{-\sqrt{\mu^2 - \omega^2}r}}{r} \quad \text{as } r \to \infty, \quad (7)$$

where $\nu = M(2\omega^2 - \mu^2)/\sqrt{\mu^2 - \omega^2}$. In the small-mass limit $M\mu \ll 1$ these solutions are well approximated by a hydrogenic spectrum [28,31] with angular separation constant $\lambda = \ell(\ell + 1)$ and frequency

$$\omega \sim \mu - \frac{\mu}{2} \left( \frac{M\mu}{\ell + n + 1} \right)^2 + \frac{i}{\gamma\ell} \left( \frac{am}{M} - 2\mu r_+ \right) (M\mu)^{4\ell+5}, \quad (8)$$

where $n = 0, 1, 2, \dots$, and $\gamma_1 = 48$ for the dominant unstable $\ell = 1$ mode.

## B. Gravitational-wave emission

For a real scalar, the condensate is a source of GWs. For a monochromatic source with frequency $\omega_R$, one can easily see by, plugging the solution (5) into the stress-energy tensor (4), that the scalar field sources GWs with frequency $2\omega_R$. In the fully relativistic regime, gravitational radiation can be computed using the Teukolsky formalism [56]. This calculation is described in detail here (see also [18,42]).

In the Teukolsky formalism, gravitational radiation is encoded in the Newman-Penrose scalar $\psi_4$, which can be decomposed as

$$\psi_4(t, r, \Omega) = \sum_{\ell m} \rho^4 \int_{-\infty}^{\infty} d\omega \sum_{\ell m} R_{\ell m \omega}(r) {}_{-2}S_{\ell m \omega}(\Omega) e^{-i\omega t}, \quad (9)$$

where $\rho = (r - ia \cos\vartheta)^{-1}$. The radial function $R(r)$ satisfies the inhomogeneous equation

$$\Delta^2 \frac{d}{dr} \left( \Delta^{-1} \frac{dR}{dr} \right) + \left( \frac{K^2 + 4i(r-M)K}{\Delta} - 8i\omega r - \lambda \right) R = T_{\ell m \omega}, \quad (10)$$

where again we omit angular indices for simplicity, $K \equiv (r^2 + a^2)\omega - am$, $\lambda \equiv A_{s\ell m} + a^2\omega^2 - 2am\omega$, and $A_{s\ell m}$ are the angular eigenvalues. The source term $T_{\ell m\omega}$ is given by

$$T_{\ell m \omega} \equiv \frac{1}{2\pi} \int d\Omega dt {}_2S_{\ell m} T e^{i\omega t}, \quad (11)$$

where $T$ is related to the scalar field stress-energy tensor (4) and can be found in [56,57].

To solve the radial equation (10) we use a Green-function approach. The Green function can be found by considering two linearly independent solutions of the homogeneous Teukolsky equation (10), with the following asymptotic behavior (see e.g. [57]):

$$R^H \to \begin{cases} \Delta^2 e^{-ikr^*} & \text{for } r \to r_+, \\ r^3 B_{\text{out}} e^{i\omega r^*} + r^{-1} B_{\text{in}} e^{-i\omega r^*} & \text{for } r \to +\infty, \end{cases} \quad (12)$$
$$R^{\infty} \rightarrow \begin{cases} A_{\text{out}} e^{ikr^{*}} + \Delta^{2} A_{\text{in}} e^{-ikr^{*}} & \text{for } r \rightarrow r_{+}, \\ r^{3} e^{i\omega r^{*}} & \text{for } r \rightarrow +\infty, \end{cases} \quad (13)$$

where $k = \omega - m\Omega_H$, $\{A, B\}_{\text{in,out}}$ are constants, and the tortoise coordinate is defined as

$$r^* = r + \frac{2Mr_+}{r_+ - r_-} \ln \frac{r - r_+}{2M} - \frac{2Mr_-}{r_+ - r_-} \ln \frac{r - r_-}{2M}. \quad (14)$$

Imposing ingoing boundary conditions at the horizon and outgoing boundary conditions at infinity, one finds that the solution of Eq. (10) is given by [57]

$$R = \frac{1}{W} \left\{ R^{\infty} \int_{r_{+}}^{r} dr' \frac{R^{H} T_{\ell m \omega}}{\Delta^2} + R^{H} \int_{r}^{\infty} dr' \frac{R^{\infty} T_{\ell m \omega}}{\Delta^2} \right\}, \quad (15)$$

where the Wronskian $W = (R^{\infty}\partial_r R^H - R^H\partial_r R^{\infty})/\Delta$ is a constant by virtue of the homogeneous Teukolsky equation (10). Using Eqs. (13) and (12) one finds

$$W = 2i\omega B_{\text{in}}. \quad (16)$$

At infinity the solutions reads

$$R(r \to \infty) \to \frac{r^3 e^{i\omega r^*}}{2i\omega B_{\text{in}}} \int_{r_+}^{\infty} dr' T_{\ell m \omega} \frac{R^H}{\Delta^2} \equiv \tilde{Z}^{\infty} r^3 e^{i\omega r^*}. \quad (17)$$

Since the frequency spectrum of the source $T_{\ell m\omega}$ is discrete with frequency $\tilde{\omega} = \pm 2\omega_{\ell mn}$ and $\tilde{m} = \pm 2m$, where $\omega_{\ell mn}$ are the scalar field eigenfrequencies, $\tilde{Z}^{\infty}$ takes the form

$$\tilde{Z}^{\infty} = \sum_{\ell \tilde{m} n} \delta(\omega - \tilde{\omega}) Z_{\ell \tilde{m} \omega}^{\infty}, \quad (18)$$

and at $r \to \infty$, $\psi_4$ is given by

$$\psi_4 = \frac{1}{r} \sum_{\ell \tilde{m} n} Z_{\ell \tilde{m} \tilde{\omega}-2}^{\infty} Y_{\ell \tilde{m} \tilde{\omega}} e^{i\tilde{\omega}(r^* - t)}. \quad (19)$$

At infinity the Newman-Penrose scalar can be written as

$$\psi_4 = \frac{1}{2}(\ddot{h}_+ - i\dot{h}_\times), \quad (20)$$

where $h_+$ and $h_\times$ are the two independent GW polarizations. The energy flux carried by these waves at infinity is given by

$$\frac{dE}{dtd\Omega} = \frac{r^2}{16\pi} (\dot{h}_+^2 + \dot{h}_\times^2). \quad (21)$$

Using Eqs. (19) and (20) we get

<figure><img src="image_3.png" /><figcaption>FIG. 2. Flux for $\ell = m = 1$ and taking the first two leading order terms in the flux $\tilde{\ell} = \tilde{m} = 2$ and $\tilde{\ell} = 3, \tilde{m} = 2$ as a function of the scalar mass and for the spin computed at the superradiant threshold (25). The numerical results computed in this work are compared with the analytic formula obtained in [18], labeled “Brito+”, and the one obtained in [17], labeled “Arvanitaki+”.</figcaption></figure>

$$\frac{dE}{dt} = \sum_{\ell \tilde{m} n} \frac{1}{4\pi\tilde{\omega}^2} |Z_{\ell \tilde{m} \tilde{\omega}}^{\infty}|^2. \quad (22)$$

We note that $|Z_{\ell \tilde{m} \tilde{\omega}}^{\infty}| \propto M_S/M^2$, where $M_S$ is the total mass of the scalar cloud:

$$M_S = \int T_i^t \sqrt{-g} dr d\theta d\varphi, \quad (23)$$

and $\sqrt{-g} = (r^2 + a^2 \cos^2 \theta) \sin \theta$ is the Kerr metric determinant. Here we neglected the energy flux at the horizon, which in general is subdominant [58]. In fact, we will only need to compute radiation at the superradiant threshold, where the flux at the horizon—being proportional to $k = (\omega - m\Omega_H)$—vanishes exactly [59].

Figure 2 shows the dominant GW energy flux computed numerically within the perturbative framework described above. Our results are compared to the analytic results of Refs. [17,18]. The flat-space approximation adopted in [17] underestimates the flux by some orders of magnitude, especially when $\mu M \ll 0.3$, for any spin. Likewise, the Schwarzschild approximation adopted in [18] overestimates the GW flux. To improve on both approximations, in the rest of this work we will use the numerical results, which are valid in the entire $(\chi, \mu M)$ plane and agree with those of [42].

### C. Evolution of the superradiant instability and of the BH-condensate system

Current nonlinear evolutions are unable to probe the development of the instability in the scalar case [48]. However, since the time scales of both the superradiant instability and the GW emission are much longer than the dynamical time scale of the BH, the evolution of the
BH-condensate system can be studied within a quasiadiabatic approximation [18]. The scalar field can be considered almost stationary, and its backreaction on the geometry neglected, as long as the scalar stress-energy tensor is small compared to the BH energy density [18].

Recent nonlinear evolutions by East and Pretorius in the spin-1 case [34,35], where the instability develops more rapidly, lend support to an adiabatic treatment of the evolution of the field. The evolution happens in two steps characterized by very different time scales. First a scalar condensate grows around the BH until the superradiant condition is saturated; then the condensate is dissipated through GW emission. Neglecting accretion for simplicity, the evolution of the system is governed by the Eq. [18]

$$ \left\{ \begin{aligned} \dot{M} &= -\dot{E}_S, \\ \dot{M} + \dot{M}_S &= -\dot{E}, \\ \dot{j} &= -m\dot{E}_S/\omega_R, \\ \dot{j} + \dot{J}_S &= -m\dot{E}/\omega_R, \end{aligned} \right. \qquad (24) $$

where $\dot{E}_S = 2M_S\omega_I$ is the scalar energy flux extracted from the horizon through superradiance. In the above equations, we have used the fact that—for a single ($\ell, m$) mode—the GW angular momentum flux is $m\dot{E}/\omega_R$ and that the angular momentum flux of the scalar field extracted at the horizon is $m\dot{E}_S/\omega_R$.

The system (24) shows that for a superradiantly unstable state ($\omega_I > 0$) the instability will cause the BH to transfer mass and spin to the scalar field until the system reaches the saturation point, given by $\omega_I = 0$, i.e., $\omega_R = m\Omega_H$.<fnref n="2" /> This process occurs on a time scale $\tau_{\text{inst}} \equiv 1/\omega_I \gg M$, and the saturation point corresponds a final BH angular momentum

$$ J_f = \frac{4mM_f^3\omega_R}{m^2 + 4M_f^2\omega_R^2} < J_i, \qquad (25) $$

where $J_{i/f}$, $M_{i/f}$ are the initial/final BH angular momentum and mass, respectively. The system (24) also shows that the variation of the BH mass $\delta M$ is related to the variation of the BH angular momentum $\delta J$ by $\delta M = \frac{\omega_R}{m}\delta J$, which implies

$$ M_f = M_i - \frac{\omega_R}{m}(J_i - J_f). \qquad (26) $$

When the instability saturates, the total mass of the scalar cloud is roughly given by $M_S^{\text{max}} \sim M_i - M_f$, namely

$$ M_S^{\text{max}} \sim \frac{J_i \omega_R}{m} - \frac{4 M_f^3 \omega_R^2}{m^2 + 4 M_f^2 \omega_R^2} \approx \frac{J_i \omega_R}{m}, \qquad (27) $$

where the last step is valid when $M_f\omega_R \ll 1$.

After the superradiant phase, the mass and the angular momentum of the BH remain constant [cf. Eq. (24)], whereas the scalar field is dissipated through the emission of GWs<fnref n="3" /> as given by Eq. (22). We neglect GW absorption at the event horizon—which is always subdominant [58]—and GW emission due to the transition of bosons between different energy levels, which is also a subdominant effect as long as the condensate is mostly populated by a single level [19]. By using again Eq. (24), after the superradiant phase we get

$$ \dot{M}_S = -\frac{dE}{dt} = -\frac{d\tilde{E} M_S^2}{dt M_f^2}, \qquad (28) $$

where we used the fact that $|Z_{\ell m\omega}^{\infty}|^2 \propto M_S^2$ to factor out the dependence on $M_S(t)$, and we defined $\frac{d\tilde{E}}{dt} \equiv \frac{dE}{dt} \frac{M_f^2}{M_S^2}$. This quantity is shown in Fig. 2 and it is constant after the superradiant phase, since it depends only on the final BH mass and spin. Therefore, setting $t = 0$ to be the time at which the superradiant phase saturates, the above equation yields

$$ M_S(t) = \frac{M_S^{\text{max}}}{1 + t/\tau_{\text{GW}}}, \qquad (29) $$

where $M_S^{\text{max}}$ is the mass of the condensate at the end of the superradiant phase [cf. Eq. (27)] and

$$ \begin{aligned} \tau_{\text{GW}} &\approx M_f \left( \frac{d\tilde{E}}{dt} \frac{M_S^{\text{max}}}{M_f} \right)^{-1} \\ &\approx 8 \times 10^5 \text{ yr} \left[ \frac{M_f}{10^6 M_{\odot}} \right] \left[ \frac{10^{-11}}{d\tilde{E}/dt} \right] \left[ \frac{0.2 M_f}{M_S^{\text{max}}} \right] \end{aligned} \qquad (30) $$

is the gravitational radiation time scale.

Finally, we note that the self-gravity of the boson cloud will cause the GW frequency to change slightly as the cloud dissipates via GWs [19,21]. The estimates given in Eq. (28) of Ref. [19] and Appendix E of Ref. [21] suggest that, for scalar fields, this small change should not affect current continuous-wave searches. Taking these estimates and the duration of the signal of Figs. 3 and 4 for resolved events, one can see that for both LIGO and LISA a vast majority of

<footnote n="2">²Fully nonlinear evolutions of a charged scalar field around a charged BH enclosed by a reflecting mirror [50,60,61] or in anti-de Sitter spacetime [51] have shown that the end state for this system indeed consists of a scalar condensate around a charged BH saturating the superradiant condition. East and Pretorius reached the same conclusion for massive spin-1 fields [34,35]. For complex fields, truly stationary metric solutions of the field equations describing a boson condensate saturating the superradiant condition around spinning BH have been explicitly shown to exist [62–64].</footnote>

<footnote n="3">³In the language of [19] this process corresponds to the “axion + axion → graviton” annihilation process. In our notation, their “occupation number” is $N = M_S/m_s$.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Gravitational radiation time scale, instability time scale, and the signal duration $\Delta t$ [defined in Eq. (33)] for detectable LISA sources and for different boson masses.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Gravitational radiation time scale, instability time scale, and the signal duration $\Delta t$ [defined in Eq. (33)] for detectable LIGO sources and for different boson masses. Dashed lines represent extragalactic sources and bold lines represent Galactic sources.</figcaption></figure>
the sources will have a small positive frequency drift $\dot{f} \ll 10^{-9}$ Hz/s, which is the current upper limit on the frequency time derivative of the latest all-sky search from LIGO [65]. However, even though this frequency drift should be very small and undetectable for most sources, the positive frequency time derivative of GWs from boson clouds could be used to distinguish them from other continuous sources, such as rotating neutron stars, which have a negative frequency drift [66].

### D. Instability and gravitational radiation time scales

As discussed above, the basic features of the evolution of the BH superradiant instability in the presence of light bosons can be understood as a two-step process, governed by two different time scales. The first time scale is the typical e-folding time of the superradiant instability given by $\tau_{\text{inst}} \equiv 1/\omega_I$, where in the $M\mu \ll 1$ limit, $\omega_I$ is the imaginary part of Eq. (8). The boson condensate grows over the time scale $\tau_{\text{inst}}$ until the superradiant condition is saturated. Subsequently, the condensate is dissipated through GW emission over a time scale $\tau_{\text{GW}}$ given by Eq. (30). In the $M\mu \ll 1$ limit, $d\tilde{E}/dt = (484 + 9\pi^2)/23040(\mu M)^{14} \approx 0.025(\mu M)^{14}$ [18,42]. Thus, using Eqs. (8), (27), (30) and reinstating physical units, the two most relevant time scales of the system are of the order

$$ \tau_{\text{inst}} \sim 10^5 \text{ yr}(M_6^8 \mu_{17}^9 \chi)^{-1}, \quad (31) $$

$$ \tau_{\text{GW}} \sim 5 \times 10^{11} \text{ yr}(M_6^{14} \mu_{17}^{15} \chi)^{-1}, \quad (32) $$

where $M_6 = M/(10^6 M_\odot)$ and $\mu_{17} = m_s/(10^{-17} \text{ eV})$ and $\chi \ll 1$.

These relations are still a reasonably good approximation when $M\mu \sim 1$ and $\chi \sim 1$. They show that there is a clear hierarchy of time scales ($\tau_{\text{GW}} \gg \tau_{\text{inst}} \gg M$), and this is important for two reasons. First of all it is crucial that $\tau_{\text{GW}} \gg \tau_{\text{inst}}$, otherwise the boson condensate would not have time to grow. Second, the time scale hierarchy justifies the use of an adiabatic approximation to describe the evolution.

Beyond the instability and gravitational radiation time scales, from the point of view of detection it is important to estimate the distribution of signal durations $\Delta t$. For LIGO we can safely neglect accretion, because accreted matter is not expected to significantly alter the birth spin of stellar-mass BHs [67]. We can also neglect the effect of mergers, since mergers affect a very small fraction of the overall population of isolated BHs [68–72], and LIGO data already suggest that multiple mergers should be unlikely [73,74]. Therefore, for LIGO we will simply assume $\Delta t = \min(\tau_{\text{GW}}, t_0)$, where $t_0 \approx 13.8$ Gyr is the age of the Universe.

For massive BHs that radiate in the LISA band, both mergers and accretion are expected to be important [75,76]. Therefore we conservatively assume that whenever an

accretion event or a merger happens the boson-condensate signal is cut short, and for LISA we define

$$ \Delta t = \left\langle \min \left( \frac{\tau_{\text{GW}}}{N_m + 1}, t_s, t_0 \right) \right\rangle, \quad (33) $$

where the signal duration $\tau_{\text{GW}}$ in the absence of mergers and accretion is given by Eq. (30), $\langle \cdots \rangle$ denotes an average weighted by the probability distribution function of the Eddington ratios, $t_s$ is the “Salpeter” accretion time scale [Eq. (51)], and $N_m$ is the average number of mergers expected in the interval $[t - \tau_{\text{GW}}/2, t + \tau_{\text{GW}}/2]$, $t$ being the cosmic time corresponding to the cosmological redshift $z$ of the GW source. Note that this definition also enforces the obvious fact that the signal cannot last longer than the age of the Universe ($\Delta t \le t_0$). We also note that the estimates of Refs. [19,21] suggest that the close passage of a stellar-mass compact object around the massive BH could affect the boson cloud when $M\mu \ll 0.1$. This part of the parameter space is mostly irrelevant for our results, and so we neglect this contribution. Moreover, estimates of the rates of extreme mass-ratio inspirals predict at most a few hundred such close passages per Gyr per galaxy [46]. Therefore, the average timescale between these events is $\gtrsim 10^7$ yr. This is comparable with the accretion timescale [Eq. (51)], which we have already taken into account. Thus, we expect our results to be robust against inclusion of this effect. In addition, stars and compact objects could, in principle, affect the boson cloud also at larger orbital distances, comparable to the peak of the cloud $R \sim 4M/(M\mu)^2$ [18]. This could also become relevant for $M\mu \ll 0.1$, but even in this case passages of stars at $R \sim 1000M$ or larger are expected to be quite rare. Indeed, tidal disruption of stars are about $10^{-5}$ per yr per galaxy [77], hence stars at distances $R \sim 1000M$ from the BH should only appear roughly every $10^5$ yr. We have checked that even if we include this effect by adding an extra time scale $\sim 10^5$ yr to Eq. (33), the background and the resolved event rates would only decrease by about an order of magnitude (and only for $M\mu \ll 0.1$), thus leaving our conclusions unchanged.

Figures 3 and 4 show histograms of $\tau_{\text{inst}}$, $\tau_{\text{GW}}$ and $\Delta t$ for resolvable sources with SNR $\rho \ge 8$ [cf. Eq. (34)]. When computing the SNR, we use an observation time $T_{\text{obs}} = 2$ yr for LIGO and $T_{\text{obs}} = 4$ yr for LISA. We adopt the LISA noise power spectral density specified in the ESA proposal for L3 mission concepts [11] and the design sensitivity of Advanced LIGO [78]. The events are binned by gravitational radiation time scale $\tau_{\text{GW}}$, instability time scale $\tau_{\text{inst}}$, and signal duration $\Delta t$, as defined in Eq. (33). For concreteness, in the plot we focus on the most optimistic astrophysical model, and we neglect the confusion noise due to the stochastic background produced by these sources (cf. [79]). For LIGO we show both Galactic and extragalactic sources.

The signal duration $\Delta t$ is typically equal to the gravitational radiation time scale $\tau_{\text{GW}}$, and (as anticipated) much longer than the instability time scale $\tau_{\text{inst}}$.
Since for LIGO we neglect the effects of mergers and accretion, the only visible difference between $\Delta t$ and $\tau_{GW}$ is due to the fact that we cut off the signal when its typical time scale is longer than the age of the Universe (i.e., as mentioned above, we set $\Delta t = t_0$ if $\tau_{GW} > t_0$). For LISA there are more subtle effects related to accretion and mergers [cf. Eq. (33)], but Figs. 3 and 4 demonstrate that the signal duration $\Delta t$ is always much longer than the instability time scale $\tau_{inst}$, as suggested by the rough estimates of Eqs. (31) and (32).

### E. Gravitational waveform

Since the GW signal from boson condensates is quasi-monochromatic, we can compute the (average) signal-to-noise ratio (SNR) as [80,81]

$$ \rho \approx \left\langle \frac{h \sqrt{T_{\text{overlap}}}}{\sqrt{S_h(f)}} \right\rangle, \quad (34) $$

where $h$ is the root-mean-square (rms) strain amplitude; $S_h(f)$ is the noise power spectral density at the (detector-frame) frequency $f$ of the signal, which is related to the source-frame frequency $f_s \equiv \omega/(2\pi)$ by $f = f_s/(1+z)$ ($z$ being the redshift); $T_{\text{overlap}}$ is the overlap time between the observation period $T_{\text{obs}}$ and the signal duration $\Delta t(1+z)$ [in the detector frame, hence the factor $1+z$ multiplying the signal duration $\Delta t$ in the source frame]; and $\langle \dots \rangle$ denotes an average over the possible overlap times. In practice, when our astrophysical models predict that a signal should overlap with the observation window, we compute this average by randomizing the signal's starting time with uniform probability distribution in the interval $[-\Delta t(1+z), T_{\text{obs}}]$ (where we assume, without loss of generality, that $t=0$ is the starting time of the observation period).

Coherent searches for almost-monochromatic sources are computationally expensive, and normally only feasible when the intrinsic parameters of the source and its sky location are known. For all-sky searches, where the properties and location of the sources are typically unknown, it is more common to use semicoherent methods, where the signal is divided in $N$ coherent segments with time length $T_{\text{coh}}$. The typical sensitivity threshold, for signals of duration $\Delta t(1+z) \gg T_{\text{obs}}$, is [cf. e.g. [66]]

$$ h_{\text{thr}} \approx \frac{25}{N^{1/4}} \sqrt{\frac{S_h(f)}{T_{\text{coh}}}}, \quad (35) $$

where $h_{\text{thr}}$ is the minimum rms strain amplitude detectable over the observation time $N \times T_{\text{coh}}$. This criterion was used, for example, in [19]. In the following we consider both cases (a full coherent search and a semicoherent method) in order to bracket uncertainties due to specific data analysis choices. For the semicoherent searches we only consider events for which $\Delta t(1+z) \gg T_{\text{obs}}$ [since the threshold given by Eq. (35) only holds for long-lived signals].

A useful quantity to compare the sensitivity of different searches independently of the data-analysis technique and the quality and amount of data is the so-called “sensitivity depth,” defined by [82]

$$ \mathcal{D}(f) = \frac{\sqrt{S_h(f)}}{h_{\text{thr}}}. \quad (36) $$

For example, the average sensitivity depth of the last Einstein@HOME search was $\mathcal{D} \approx 35 \text{ Hz}^{-1/2}$ [83].

To compute $h$, we first use Eqs. (9), (19) and (20) to get a combination of the two GW polarizations,

$$ H \equiv h_+ - ih_\times = -\frac{2}{\tilde{\omega}^2} \sum_{\ell \tilde{m} n} Z_{\ell \tilde{m} \tilde{\omega}-2}^{\infty} Y_{\ell \tilde{m} \tilde{\omega}} e^{i\tilde{\omega}(r^* - t)}. \quad (37) $$

In the following we will omit the sum over $\ell \tilde{m} n$ for ease of notation. Let us focus on a single scalar field mode.<fnref n="4" /> If the scalar field has azimuthal number $m$ and real frequency $\omega_R$, the GW emitted by the scalar cloud will have azimuthal number $\tilde{m} = \pm 2m$ and frequency $\tilde{\omega} = \pm 2\omega_R$. Defining $Z^{\infty} = |Z|e^{-i\phi}$, where $|Z|$ and $\phi$ are both real, we have

$$ H = -\frac{2|Z|}{\tilde{\omega}^2 r} \left( -2Y_{\ell \tilde{m} \tilde{\omega}} e^{i[\tilde{\omega}(r^* - t) + \phi]} + -2Y_{\ell - \tilde{m} - \tilde{\omega}} e^{-i[\tilde{\omega}(r^* - t) + \phi]} \right), \quad (38) $$

where we used the fact that $Z_{\ell-\tilde{m}-\tilde{\omega}}^{\infty} = Z_{\ell \tilde{m} \tilde{\omega}}^{\infty}$. Since $Y_{\ell \tilde{m} \tilde{\omega}}(\theta, \phi) = S_{\ell \tilde{m} \tilde{\omega}}(\theta)e^{im\phi}$ and $S$ is a real function for real $\tilde{\omega}$, we get

$$ h_+ = \Re(H) = -\frac{2|Z|}{\tilde{\omega}^2 r} (-2S_{\ell \tilde{m} \tilde{\omega}} + -2S_{\ell-\tilde{m}-\tilde{\omega}}) \times \cos[\tilde{\omega}(r^* - t) + \phi + \tilde{m}\phi], \quad (39) $$

$$ h_\times = \Im(H) = -\frac{2|Z|}{\tilde{\omega}^2 r} (-2S_{\ell \tilde{m} \tilde{\omega}} - -2S_{\ell-\tilde{m}-\tilde{\omega}}) \times \sin[\tilde{\omega}(r^* - t) + \phi + \tilde{m}\phi]. \quad (40) $$

The GW strain measured at the detector is

$$ h = h_+ F_+ + h_\times F_\times, \quad (41) $$

where $F_{+,\times}$ are pattern functions that depend on the orientation of the detector and the direction of the source. To get the rms strain of the signal we angle-average over source and detector directions and use $\langle F_+^2 \rangle = \langle F_-^2 \rangle = 1/5$, $\langle F_+ F_- \rangle = 0$, $\langle |S_{\ell \tilde{m} \tilde{\omega}}|^2 \rangle = 1/(4\pi)$ and $\langle \cos^2[\tilde{\omega}(r^* - t) + \phi + \tilde{m}\phi] \rangle = \langle \sin^2[\tilde{\omega}(r^* - t) + \phi + \tilde{m}\phi] \rangle = 1/2$. We then obtain

<footnote n="4">⁴In this work we will focus on the mode with the smallest instability time scale $\ell = m = 1$, which should be the dominant source of GW radiation [19].</footnote>
$$h \approx \langle h^2 \rangle^{1/2} = \left( \frac{2|Z|^2}{5\pi\tilde{\omega}^4 r^2} \right)^{1/2} = \left( \frac{4\dot{E}}{5\tilde{\omega}^2 r^2} \right)^{1/2}, \quad (42)$$

where $\dot{E}$ is given in Eq. (22), which for a single scalar mode reads $\dot{E} = \sum_{\ell} |Z_{\ell}|^2/(2\pi\tilde{\omega}^2)$. Finally, let us factor out the BH mass and the mass of the scalar condensate: $|Z| = A(\chi, \mu M)(M\tilde{\omega})^2 M_S/M^2$, where $A(\chi, \mu M)$ is a dimensionless quantity. The final expression for the rms strain reads

$$h = \sqrt{\frac{2}{5\pi}} \frac{M}{r} \frac{M_S}{M} A(\chi, \mu M). \quad (43)$$

We conservatively assume that the GWs observed at the detector are entirely produced after the saturation phase of the instability. Therefore, we compute $h$ using the final BH mass and spin, as computed in Eqs. (26) and (25), respectively. Larger initial spins imply that a larger fraction of the BH mass is transferred to the scalar condensate [cf. Eq. (27)]. So, for a given scalar field mass and initial BH mass, the strain grows with the initial spin.

Equation (43) is valid for any interferometric detector for which the arms form a 90-degree angle, such as Advanced LIGO. For a triangular LISA-like detector the arms form a 60-degree angle, and we must multiply all amplitudes by a geometrical correction factor $\sqrt{3}/2$ [5,84]. Additionally, since we sky-average the signal, we will use an effective non-sky-averaged noise power spectral density, obtained by multiplying LISA's sky-averaged $S_h$ by 3/20 [85]. The analysis presented below takes into account these corrective factors.

### F. Cosmological effects

Since some sources can be located at non-negligible redshifts, the root-mean-square strain amplitude of Eqs. (42) and (43) must be corrected to take into account cosmological effects, which affect the propagation of the waves to the detector [86]. These effects have two main consequences.

First, the frequency $f$ of the signal as measured at the detector's location ("detector frame") is redshifted with respect to the emission frequency $f_s$ in the "source-frame", i.e., $f = f_s/(1+z)$.

Second, in the strain amplitude given by Eq. (43), the distance $r$ to the detector should be interpreted as the comoving distance, which for a flat Friedmann-Lemaitre-Robertson-Walker model is given by

$$D_c(z) = D_H \int_0^z \frac{dz'}{\sqrt{\Delta(z')}}, \quad (44)$$

where $\Delta(z) = \Omega_M(1+z)^3 + \Omega_\Lambda$, $D_H$ is the Hubble distance, $\Omega_M$ is the dimensionless matter density and $\Omega_\Lambda$ is the dimensionless cosmological constant density. All other quantities (masses, lengths and frequencies) in Eq. (43) should be instead be interpreted as measured by an observer in the source frame.

Alternatively, one might wish to use quantities measured by an observer at the detector's location to compute the strain amplitude of Eq. (43). Detector-frame quantities are related to source-frame ones by powers of $(1+z)$, namely all quantities with dimensions $[\text{mass}]^p$ (in our geometrized units $G = c = 1$) are multiplied by the factor $(1+z)^p$, e.g. masses are multiplied by $(1+z)$ ("redshifted masses"), frequencies are divided by the same factor ("redshifted frequencies"), while the comoving distance is multiplied by a factor $(1+z)$, thus becoming the luminosity distance $D_L = D_c(1+z)$. Since the strain amplitude of Eq. (43) is dimensionless, that equation yields the same result when using detector-frame quantities as when using source-frame ones.

The typical distance up to which BH-condensate sources are detectable can be estimated by defining an "angle-averaged range" $D_{\text{range}}$ as the luminosity distance at which either the SNR $\rho(D_{\text{range}}) = 8$ [cf. Eq. (34)] for coherent searches, or $h(D_{\text{range}}) = h_{\text{thr}}$ for semicoherent searches [cf. Eq. (35)].

In Fig. 5 we show $D_{\text{range}}$ for both LISA and LIGO at design sensitivity under different assumptions on the initial BH spin. The left panels refer to single coherent observation with $T_{\text{obs}} = 4$ yr for LISA ($T_{\text{obs}} = 2$ yr for Advanced LIGO), whereas the right panels refer to a (presumably more realistic) semicoherent search with $N = 121$ coherent segments of duration $T_{\text{coh}} = 250$ hr. In the more optimistic case, sources are detectable up to cosmological distances of ~20 Gpc (~2 Gpc) if the BH is nearly extremal and the boson mass is in the optimal mass range $m_s \sim 10^{-17}$ eV ($m_s \sim 10^{-13}$ eV) for LISA (LIGO). For the semicoherent search, $D_{\text{range}}$ is reduced by roughly one order of magnitude, with a maximum detector reach ~2 Gpc and ~200 Mpc for LISA and Advanced LIGO, respectively.

### III. MASSIVE BLACK HOLE POPULATION MODELS

An assessment of the detectability of GWs from super-radiant instabilities requires astrophysical models for the massive BH population. In this section we describe the models adopted in our study, and in particular our assumptions on (A) the mass and spin distribution of isolated massive BHs, (B) their Eddington ratio distribution, and (C) their merger history.

#### A. Mass and spin distribution of isolated black holes

Let $n$ be the comoving-volume number density of BHs. For the mass and spin distribution of isolated BHs we consider:

(A.1) A model where $d^2n/(d \log_{10} M d\chi)$ is computed using the semianalytic galaxy formation model of [43] (with later improvements described in [76,87,88]). This distribution is redshift-dependent
<figure><img src="image_6.png" /><figcaption>FIG. 5. Angle-averaged range $D_{\text{range}}$ for LISA (top) and Advanced LIGO at design sensitivity (bottom) computed for selected initial BH spin ($\chi_i = 0.998, 0.95, 0.7$). Left panels: the range is computed using a coherent search over an observation time $T_{\text{obs}} = 4$ yr (for LISA) and $T_{\text{obs}} = 2$ yr (for LIGO). Right panels: we assume a semicoherent search with $N = 121$ coherent segments of duration $T_{\text{coh}} = 250$ hr.</figcaption></figure>

and skewed toward large spins, at least at low masses (cf. [76]). It also has a negative slope $dn/d\log_{10}M \propto M^{-0.3}$ for BH masses $M < 10^7 M_\odot$, which is compatible with observations (cf. [76], Fig. 7). The normalization is calibrated so as to reproduce the observed $M-\sigma$ and $M-M_*$ scaling relations of [89], where $\sigma$ is the galaxy velocity dispersion and $M_*$ is the stellar mass. We also account for the bias due to the resolvability of the BH sphere of influence [90,91]. Because of the slope, normalization and spin distribution, this model is *optimistic*.

### (A.2) An analytic mass function [46,47]

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

which we use for redshifts and BH masses in the range $10^4 M_\odot < M < 10^7 M_\odot$ and $z < 3$. For $M > 10^7 M_\odot$ we use a mass distribution with normalization 10 times lower than the optimistic one. For this model we use a *uniform* distribution of the initial spins $\chi \in [0, 1]$. Because of the lower normalization and the spin distribution, this model is *less optimistic*.

### (A.3) An analytic mass function

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

which we use again for $10^4 M_\odot < M < 10^7 M_\odot$ and $z < 3$, whereas for $M > 10^7 M_\odot$ we use a mass distribution with normalization 100 times lower than the optimistic one. For this model we also consider a *uniform* distribution of the initial spins $\chi \in [0, 1]$. Because of the normalization, slope and spin distribution, this model is *pessimistic*.

## B. Black hole mergers

Our standard choice for BH mergers is to compute the comoving-volume number density $n_m$ of mergers per (logarithmic) unit of total mass $M_{\text{tot}} = M_1 + M_2$, unit redshift and (logarithmic) unit of mass ratio $q = M_2/M_1 \le 1$, i.e.,

$$ \nu(M_{\text{tot}}, z, q) = \frac{d^3 n_{\text{m}}}{d\log_{10}M_{\text{tot}}dzd\log_{10}q}, \quad (47) $$

from the semianalytic model of [43].
We can then estimate the average number of mergers (between z and z + dz) for a BH of mass M as

$$dN_m(M, z) = \frac{\mu(M, z)}{\phi(M, z)} dz. \quad (48)$$

Here

$$\phi(M, z) \equiv \frac{dn}{d\log_{10} M} = \int \frac{d^2 n}{d\log_{10} M d\chi} d\chi \quad (49)$$

is the (isolated BH) mass function, and

$$\mu(M_{\text{tot}}, z) \equiv \frac{d^2 n_{\text{merger}}}{d \log_{10} M_{\text{tot}} dz} = \int_{q>q_c} \nu d \log_{10} q,$$

where $q_c$ is the critical mass ratio above which we assume mergers make an impact. In practice, most BH mergers in our semianalytic models have $q \gtrsim 0.01-0.001$ (especially in the LISA band, cf. [92]), so our results are robust against the exact choice of $q_c$. Nevertheless, to be on the conservative side, we set $q_c = 0$. A larger $q_c$ would produce a slightly lower BH merger number and, in turn, a slightly higher number of boson-condensate sources, under the conservative assumption that mergers destroy the boson cloud. We can then compute the average number of mergers experienced by a BH of mass M in the redshift interval $[z_1, z_2]$ as

$$N_m = \int_{z_1}^{z_2} \frac{dN_m}{dz} dz. \quad (50)$$

Note that the number of mergers depends on the seeding mechanisms of the massive BH population, as well as on the “delays” between the mergers of galaxies and the mergers of the BHs they host [cf. e.g. [39]].

When computing the average number of mergers $N_m$ to be used to estimate the number of boson-condensate GW events from isolated BHs, i.e., when evaluating the number of resolved events [Eq. (62) below] and the amplitude of the stochastic background [Eq. (64) below], we consider the “popIII” model of [39] (a light-seed scenario with delays). Choosing a different seed model would not alter our conclusions. However, when considering the constraints that can be placed on the boson mass by direct observations of BH coalescences by LISA, we consider all three models presented in [39] (“popIII”, “Q3” and “Q3nod”). These models correspond respectively to light seeds with delays between a galaxy merger and the corresponding binary BH merger; heavy seeds with delays; and heavy seeds with no delays; and they are chosen to bracket the theoretical uncertainties on the astrophysics of BH seed formation and BH delays.

### C. Accretion

Clearly, accretion is competitive with the superradiant extraction of angular momentum from the BH [18], so it is important to quantify its effect. We estimate the accretion time scale via the Salpeter time,

$$t_S = 4.5 \times 10^8 \text{ yr} \frac{\eta}{f_{\text{Edd}}(1-\eta)}, \quad (51)$$

where $f_{\text{Edd}}$ is the Eddington ratio for mass accretion, and the thin-disk radiative efficiency $\eta$ is a function of the spin related to the specific energy $E_{\text{ISCO}}$ at the innermost stable circular orbit [93]:

$$\eta = 1 - E_{\text{ISCO}}, \quad (52)$$

$$E_{\text{ISCO}} = \sqrt{1 - \frac{2}{3r_{\text{ISCO}}}}, \quad (53)$$

$$r_{\text{ISCO}} = 3 + Z_2 - \frac{\chi}{|\chi|} \sqrt{(3 - Z_1)(3 + Z_1 + 2Z_2)}, \quad (54)$$

$$Z_1 = 1 + (1 - \chi^2)^{1/3} [(1 + \chi)^{1/3} + (1 - \chi)^{1/3}], \quad (55)$$

$$Z_2 = \sqrt{3\chi^2 + Z_1^2}. \quad (56)$$

For the Eddington ratio $f_{\text{Edd}}$ we consider three models:

(C.1) We use the results of our semianalytic model to construct probability distribution functions for $f_{\text{Edd}}$ at different redshifts and BH masses.

(C.2) We adopt a simple model in which $f_{\text{Edd}} = 1$ for 10% of the massive BHs, and $f_{\text{Edd}} = 0$ for the remaining ones. (The choice of 10% is a reasonable estimate for the duty cycle of active galactic nuclei [94,95]).

(C.3) Finally, we consider a very pessimistic model in which all BHs have $f_{\text{Edd}} = 1$. Although unrealistic, this models maximizes the effects of accretion, and therefore it yields the most conservative lower bound for the superradiant instability time scale.

## IV. STELLAR MASS BLACK HOLE POPULATION MODELS

We now turn to a description of stellar-mass BHs, which are of interest for LIGO. Here we have to model (A) extragalactic BHs, which turn out to dominate the stochastic background of GWs from ultralight bosons, and (B) Galactic BHs, which (as pointed out in [19,20]) are dominant in terms of resolvable signals.

### A. Extragalactic BHs

In the standard scenario, stellar-mass BHs are the end products of the evolution of massive ($M \gtrsim 20 M_\odot$) stars. They form either via direct collapse of the star or via a supernova explosion followed by fallback of matter (failed supernova). This process depends on various parameters, such as stellar metallicity, rotation and interactions with a companion if the star belongs to a binary system [96–99]. In particular, the metallicity of the star determines the strength of stellar winds and can thus have a significant impact on the mass of the stellar core prior to collapse [100,101]. In
addition, BHs can grow hierarchically through multiple mergers that occur in dense stellar clusters [73,74,102, 103]. This process is expected to leave an imprint on the distribution in the mass-spin plane: while BHs grow in mass via mergers their spins converge to values around ~0.7 with little or no support below ~0.5 [73–75].

In this work we consider only BH formation from core collapse of massive stars. We use the analytic fits for the BH mass as a function of initial stellar mass and metallicity from [104], embedded in the semianalytic galaxy evolution model from [105]. In particular, the latter model describes the production of metals by stars [106] and the evolution of the metallicity of the interstellar medium, which is inherited by the stars that form there. The extragalactic BH formation rate as a function of mass and redshift reads

$$ \frac{d\dot{n}_{\text{eg}}}{dM} = \int d\mathcal{M}_{\star} \psi[t - \tau(\mathcal{M}_{\star})] \phi(\mathcal{M}_{\star}) \delta[\mathcal{M}_{\star} - g^{-1}(M)], \quad (57) $$

where $\tau(\mathcal{M}_{\star})$ is the lifetime of a star with mass $\mathcal{M}_{\star}$, $\phi(\mathcal{M}_{\star})$ is the stellar initial mass function, $\psi(t)$ denotes the cosmic star formation rate (SFR) density and $\delta$ is the Dirac delta. We use the fit to the cosmic SFR described in [107], calibrated to observations [108,109]. We adopt a Salpeter initial mass function $\phi(\mathcal{M}_{\star}) \propto \mathcal{M}_{\star}^{-2.35}$ [110] in the mass range $\mathcal{M}_{\star} \in [0.1–100] M_{\odot}$ and use the stellar lifetimes from [111]. The initial stellar mass $\mathcal{M}_{\star}$ and BH mass $M$ are related by the function $M = g(\mathcal{M}_{\star})$, which can be (implicitly) redshift-dependent (through its dependence on stellar metallicity), and which we take from the “delayed” model of [104].

## B. Galactic BHs

Resolvable signals are expected to be dominated by Galactic stellar-mass BHs [19]. We estimate the present-day mass function of these BHs as

$$ \frac{dN_{\text{MW}}}{dM} = \int dt \frac{\text{SFR}(z)}{\mathcal{M}_{\star}} \frac{dp}{d\mathcal{M}_{\star}} \left| \frac{dM}{d\mathcal{M}_{\star}} \right|^{-1}, \quad (58) $$

where $N_{\text{MW}}$ denotes the number of BHs in the Galaxy, $dp/d\mathcal{M}_{\star}$ is the normalized Salpeter initial mass function (i.e., the probability of forming a star with mass between $\mathcal{M}_{\star}$ and $\mathcal{M}_{\star} + d\mathcal{M}_{\star}$), and SFR($z$) denotes the SFR of Milky-Way type galaxies as a function of $z$ [109,112]. The integration is over all cosmic times till the present epoch. The (differential) relation between BH mass and initial stellar mass $dM/d\mathcal{M}_{\star}$ is taken from the “delayed” model of [104], and is also a function of redshift via the metallicity. For the latter, we use the results of [113] to describe its evolution with cosmic time. We then “spread” $dN_{\text{MW}}/dM$ throughout the Galaxy in order to obtain a (differential) density $dn_{\text{MW}}/dM$, by assuming that the latter is everywhere proportional to the (present) stellar density. To this purpose, we describe the Galaxy by a bulge + disk model, where the bulge follows a Hernquist profile [114]

with mass ~2 × 10¹⁰ $M_{\odot}$ and scale radius ~1 kpc [115], and the disk is described by an exponential profile with mass ~6 × 10¹⁰ $M_{\odot}$ and scale radius ~2 kpc [116].

Since these models (for both Galactic and extragalactic BHs) do not predict the initial BH spins, we assume a uniform distribution and explore different ranges (from optimistic to pessimistic): $\chi \in [0.8, 1]$, [0.5, 1], [0, 1] and [0, 0.5].

## V. EVENT RATES FOR LISA AND LIGO

Having in hand the calculation of the GW signal of Sec. II and the astrophysical models of Secs. III and IV, we can now compute event rates for LISA and LIGO. We consider two separate classes of sources: (A) boson-condensate GW events which are loud enough to be individually resolvable, and (B) the stochastic background of unresolvable sources.

### A. Resolvable sources

In the limit in which the (detector-frame) signal duration $\Delta t(1+z)$ is small compared to the observation time $T_{\text{obs}}$, $\Delta t(1+z) \ll T_{\text{obs}}$, the number of resolvable events is proportional to the observation time [117]:

$$ N = T_{\text{obs}} \int_{\rho>8} \frac{d^2 \dot{n}}{dM d\chi} \frac{dt}{dz} 4\pi D_c^2 dz dM d\chi, \quad (59) $$

where

$$ \frac{dt}{dz} = \frac{1}{H_0 \sqrt{\Delta}(1+z)} \quad (60) $$

is the derivative of the lookback time with respect to redshift.

For long-lived sources with detector-frame duration $\Delta t(1+z) \gg T_{\text{obs}}$, the number of detections does not scale with the observation time, but rather with the “duty cyle” $\Delta t/t_f$, where $t_f \equiv n/\dot{n}$ is the formation time scale of the boson condensate. For example, if BHs form a boson condensate only once in their cosmic history, $t_f$ is the age of the Universe $t_0 \approx 13.8$ Gyr. This duty cycle has the same meaning as the duty cycle of active galactic nuclei: it accounts for the fact that, at any given time, only a fraction of the BH population will be emitting GWs via boson condensates. Because of the ergodic theorem, this fraction is given by the average time fraction during which a BH emits GWs via boson condensates. This average time fraction is indeed the duty cycle $\Delta t/t_f$. Therefore, the number of resolved sources when $\Delta t(1+z) \gg T_{\text{obs}}$ is simply

$$ \begin{aligned} N &= \int_{\rho>8} \frac{d^2 n}{dM d\chi} \frac{\Delta t}{t_f} \frac{dV_c}{dz} dz dM d\chi \\ &= \int_{\rho>8} \frac{d^2 \dot{n}}{dM d\chi} \Delta t \frac{dV_c}{dz} dz dM d\chi, \end{aligned} \quad (61) $$

where $dV_c = 4\pi D_c^2 dD_c$.
Equations (59) and (61) can be merged into a single expression that remains valid also in the intermediate regime $\Delta t(1+z) \sim T_{\text{obs}}$. Indeed, the probability that a signal lasting a time span $\Delta t(1+z)$ (in the detector frame) overlaps with an observation of duration $T_{\text{obs}}$ is simply proportional to the sum of the two durations, $\Delta t(1+z) + T_{\text{obs}}$. This can be understood in simple geometric terms: for the signal to overlap with the observation window (which we define, without loss of generality, to extend from $t=0$ to $t=T_{\text{obs}}$), the signal's starting time should fall between $t = -\Delta t(1+z)$ and $t = T_{\text{obs}}$, i.e., in a time interval of length $\Delta t(1+z) + T_{\text{obs}}$. Therefore, we can estimate the number of observable GW events as

$$N = \int_{\rho>8} \frac{d^2\dot{n}}{dM d\chi} \left( \frac{T_{\text{obs}}}{1+z} + \Delta t \right) \frac{dV_c}{dz} dz dM d\chi. \quad (62)$$

Since $dD_c/dz = (1+z)dt/dz$, it can be easily checked this equation reduces to Eqs. (59) and (61) in the limits $\Delta t(1+z) \ll T_{\text{obs}}$ and $\Delta t(1+z) \gg T_{\text{obs}}$, respectively.

For extragalactic LIGO sources we compute $d^2\dot{n}/dM d\chi$ from the astrophysical models of Sec. IVA, while for LISA and galactic LIGO sources we compute $d^2n/dM d\chi$ as described in Secs. III and IV B and then assume $d^2\dot{n}/dM d\chi = (d^2n/dM d\chi)/t_0$. This corresponds to assuming that the boson-condensate formation time $t_f = t_0$ equals the age of the Universe, or that BHs radiate via boson condensates only once in their lifetime. This conservative assumption does not affect our results very significantly. Once a BH-boson system radiates, its spin decreases to low values, while the mass remains almost unchanged. For the BH to emit again via boson condensates, its spin must grow again under the effect of accretion or mergers. In this process, however, the BH mass also grows rapidly: for example, the simple classic estimates by Bardeen [118] imply that when a BH spins up from $\chi=0$ to $\chi=1$ via accretion, its mass increases by a factor $\sqrt{6}$. So even if new boson clouds form due to the instability of higher-$m$ modes, the instability time scales will be much larger [cf. Eq. (8)] and the GW flux will be highly suppressed [cf. Ref. [42]].

Our main results for resolvable rates are summarized in Fig. 6, Fig. 7, Table I and Table II.

In Fig. 6 we focus on optimistic models and we show how the number of individually resolvable events depends on the observation time and on the chosen data-analysis method. More specifically, for LISA we use the BH mass-spin distribution model (A.1) and accretion model (C.1), while for LIGO we adopt the optimistic spin distribution $\chi_i \in [0.8, 1]$. We bracket uncertainties around the nominal LISA mission duration of $T_{\text{obs}} = 4$ yr [11] by considering single observations with duration $T_{\text{obs}} = (2, 4, 10)$ yr. We also show rates for a (presumably more realistic) semicoherent search with 121 segment of $T_{\text{coh}} = 250$ hours coherent

<figure><img src="image_7.png" /><figcaption>FIG. 6. Number of resolved LIGO and LISA events for our optimistic BH population models as a function of the boson mass with different observation times $T_{\text{obs}}$, using both full and semi-coherent searches. Thick (thin) lines were computed with (without) the confusion noise from the stochastic background.</figcaption></figure>

integration time<fnref n="5" /> For Advanced LIGO at design sensitivity, we similarly consider single observations lasting either $T_{\text{obs}} = 2$ yr or $T_{\text{obs}} = 4$ yr, as well as a semi-coherent search with 121 segments of $T_{\text{coh}} = 250$ hours coherent integration time.

Figure 6 (together with Fig. 3 in [79]) shows that the number of resolvable events is strongly dependent on the boson mass and on the astrophysical model.

For LISA, our astrophysical populations contain mostly BHs in the mass range $10^4 M_{\odot} < M < 10^8 M_{\odot}$, and the sensitivity curve peaks around a frequency corresponding to $m_s \sim 10^{-17}$ eV [cf. Fig. 1 of [79]]. These considerations—together with the condition for having an efficient superradiant instability (namely, $M\mu \sim 0.4$ at large spin)—translate into the range $3 \times 10^{-18}$ eV $\lesssim m_s \lesssim 5 \times 10^{-17}$ eV for the mass of detectable bosonic particles in a semicoherent search.

For LIGO, our models predict that most BHs will be in the mass range $3 M_{\odot} < M < 50 M_{\odot}$, and the most sensitive frequency band corresponds to $m_s \sim 3 \times 10^{-13}$ eV [cf. Fig. 1 of [79]], translating into the range $2 \times 10^{-13}$ eV $\lesssim m_s \lesssim 3 \times 10^{-12}$ eV for the mass of detectable bosonic particles.

In order to quantify the “self-confusion” noise due to the stochastic background produced by BH-boson systems, in Fig. 6 we also display the number of resolved events that we would obtain if we omitted the confusion noise from the stochastic background (cf. Fig. 1 of [79] and Sec. V B). Neglecting the confusion noise would overestimate the number of resolvable events in LISA by one or two orders of magnitude.

<footnote n="5">⁵The number of resolved events for other choices of number of segments and coherent integration time can be obtained from Fig. 7 and expressing the sensitivity depth as $D \approx T_{\text{coh}}^{1/2} N^{1/4} 25^{-1}$ [cf. Eqs. (35) and (36)].</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Left: Number of events as a function of the sensitivity depth $\mathcal{D}$ [Eq. (36)] for selected boson masses in the LISA band and accretion model (C.1). The bottom (top) of each shadowed region correspond to the pessimistic (optimistic) model. Right: Same, but for boson masses in the LIGO band. Here the bottom (top) of each shadowed region correspond to pessimistic (optimistic) spin distributions.</figcaption></figure>

The rates computed in Figure 6 refer to our optimistic astrophysical models. As shown in [79], resolvable event rates in the most pessimistic models are about one order of magnitude lower. Nevertheless, it is remarkable that even in the most pessimistic scenario for direct detection (i.e., unfavorable BH mass-spin distributions and semicoherent search method for the signal), bosonic particles with $m_s \sim 10^{-17}$ eV ($m_s \sim 10^{-12}$ eV) would still produce around 5 (15) direct LISA (LIGO) detections of boson-condensate GW events.

In Fig. 7 we show how the number of events grows with the sensitivity depth of the search [82], as defined in Eq. (36). For LISA the number of events grows roughly with $\mathcal{D}^3$, corresponding to $T_{\text{obs}}^{3/2}$. This is expected from the fact that the number of events for sources at $\gtrsim 30$ Mpc should grow with the sensitive volume, and thus decrease with $\rho_{\text{crit}}^{-3}$, where $\rho_{\text{crit}}$ is the critical SNR for detection [119].

On the other hand, LIGO will be mostly sensitive to signals within the Galaxy. For a given boson mass and distance, $\tau_{\text{GW}} \sim h^{-2}$ and $M \sim h^{1/8}$ [cf. Eqs. (30) and (42)]. Since the Galactic stellar BH population obtained from Eq. (58) is well fitted by $dN/M \sim e^{-0.2M}$, for a fixed volume the integral in Eq. (61) goes as

$$N \sim \int_{h>h_{\text{thr}}} h^{-23/8} e^{-0.2h^{1/8}} dh \sim h_{\text{thr}}^{-15/8}, \quad (63)$$

where in the last step we took the leading order of the integral for small $h_{\text{thr}}$. From Eq. (34) one has $h_{\text{thr}} \propto T_{\text{obs}}^{-1/2}$ and therefore $N \propto T_{\text{obs}}^{15/16}$. This is in agreement with the scaling that we find.

Assuming the sensitivity depth of the last EINSTEIN@ HOME search $\mathcal{D} \approx 35 \text{ Hz}^{-1/2}$ [83] and an optimal boson mass around $m_s \sim 10^{-12.5}$ eV, we find that O1 should have detected 5 resolvable events for the optimistic spin distribution $\chi \in [0.8, 1]$, and 2 events for a uniform spin distribution $\chi \in [0, 1]$.

As pointed out in [79], these optimal boson masses may already be ruled out by upper limits from existing stochastic background searches [79]. On the other hand, the pessimistic spin distribution $\chi \in [0, 0.5]$ is still consistent with (the lack of) observations of resolvable BH-boson GW events in O1, though marginally ruled out by the O1 stochastic background upper limits [79].

<table><caption>TABLE I. Number of resolvable events in the LISA band computed including the “self-confusion” noise from the stochastic background of BH-boson condensates for different accretion models. The lower and upper bounds correspond to the pessimistic and optimistic massive BH population models, respectively. For the semicoherent search we use 121 segments of $T_{\text{coh}} = 250$ hours coherent integration time. For the coherent search, we adopt the nominal mission duration of $T_{\text{obs}} = 4$ years.</caption><thead><tr><th>m<sub>s</sub> [eV]</th><th>Search method</th><th>Accretion model</th><th>Events</th></tr></thead><tbody><tr><td rowspan="6">10<sup>-16</sup></td><td>Coherent</td><td>(C.1)</td><td>75-0</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr><tr><td>Coherent</td><td>(C.2)</td><td>75-0</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr><tr><td>Coherent</td><td>(C.3)</td><td>75-0</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr><tr><td rowspan="6">10<sup>-17</sup></td><td>Coherent</td><td>(C.1)</td><td>1329-1022</td></tr><tr><td>Semicoherent</td><td></td><td>39-5</td></tr><tr><td>Coherent</td><td>(C.2)</td><td>3865-1277</td></tr><tr><td>Semicoherent</td><td></td><td>36-4</td></tr><tr><td>Coherent</td><td>(C.3)</td><td>5629-1429</td></tr><tr><td>Semicoherent</td><td></td><td>39-5</td></tr><tr><td rowspan="6">10<sup>-18</sup></td><td>Coherent</td><td>(C.1)</td><td>17-1</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr><tr><td>Coherent</td><td>(C.2)</td><td>18-1</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr><tr><td>Coherent</td><td>(C.3)</td><td>20-0</td></tr><tr><td>Semicoherent</td><td></td><td>0</td></tr></tbody></table>
<table><caption>TABLE II. Number of resolvable events for Advanced LIGO at design sensitivity. For the semicoherent search we use 121 segments of $T_{\text{coh}} = 250$ hours coherent integration time. For the coherent search, we set $T_{\text{obs}} = 2$ years. The lower and upper bounds correspond to the pessimistic ($\chi \in [0, 0.5]$) and optimistic ($\chi \in [0.8, 1]$) spin distributions, respectively.</caption><thead><tr><th>$m_s$[eV]</th><th>Search method</th><th>Events</th></tr></thead><tbody><tr><td rowspan="2">$10^{-11.5}$</td><td>Coherent</td><td>21–2</td></tr><tr><td>Semicoherent</td><td>1–0</td></tr><tr><td rowspan="2">$10^{-12}$</td><td>Coherent</td><td>1837–193</td></tr><tr><td>Semicoherent</td><td>50–2</td></tr><tr><td rowspan="2">$10^{-12.5}$</td><td>Coherent</td><td>12556–1429</td></tr><tr><td>Semicoherent</td><td>205–15</td></tr></tbody></table>

Our results for resolvable event rates using different search techniques, mass/spin and accretion models are summarized in Tables I and II. For LISA we included “self-confusion” noise in our rate estimates, and using different accretion models does not significantly affect our results. Interestingly, even though the accretion models (C.2) and (C.3) are more pessimistic than model (C.1), they predict a slightly larger number of resolvable events for boson masses in the optimal range around $10^{-17}$ eV. This is because the self-confusion noise is lower for models (C.2) and (C.3) [cf. Sec. V B], and thus the loss in signal is more than compensated by the lower total (instrumental and self-confusion) noise floor.

## B. Stochastic background

In addition to individually resolvable sources, a population of massive BH-boson condensates at cosmological distances can build up a detectable stochastic background. This possibility is potentially very interesting, given the spread in BH masses (and, hence, in boson masses that would yield an instability) characterizing the BH population at different redshifts, but to the best of our knowledge it has not been explored in the existing literature.

The stochastic background can be computed from the formation rate density per comoving volume $\dot{n}$ as [120]

$$ \Omega_{\text{gw}}(f) = \frac{f}{\rho_c} \int_{\rho<8} d\chi dMdz \frac{dt}{dz} \frac{d^2\dot{n}}{dMd\chi} \frac{dE_s}{df_s}, \quad (64) $$

where $\rho_c = 3H_0^2/(8\pi G)$ is the critical density of the Universe, $dE_s/df_s$ is the energy spectrum in the source frame, and $f$ is the detector-frame frequency. Note that the integral is only over unresolved sources with $\rho < 8$.

For extragalactic stellar mass BHs (which are sources for LIGO), we calculate $d^2\dot{n}/dMd\chi$ based on the model of Sec. IV, while for LISA sources we use the model of Sec. III to obtain $d^2n/dMd\chi$, and then (as we did for the resolved sources) we assume $d^2\dot{n}/dMd\chi = (1/t_0)(d^2n/dMd\chi)$. As

before, this corresponds to the conservative assumption that formation of boson condensates occurs only once in the cosmic history of each massive BH.

We compute the energy spectrum as

$$ \frac{dE_s}{df_s} \approx E_{\text{GW}}\delta(f(1+z) - f_s), \quad (65) $$

where we recall that $f_s$ is the frequency of the signal in the source frame, $E_{\text{GW}}$ is the total energy radiated by the boson cloud in GWs during the signal duration $\Delta t$, and the Dirac delta is “spread out” over a frequency window of size $\sim \max[1/(\Delta t(1+z)), 1/T_{\text{obs}}]$ to account for the finite signal duration and the finite frequency resolution of the detector. As in the calculation of the rates of resolved sources, $\Delta t = \min(\tau_{\text{GW}}, t_0)$ [cf. Eq. (30)] for LIGO sources, while we account for mergers and accretion through Eq. (33) for LISA sources. Moreover, since our calculations rely on the implicit assumption that the instability reaches saturation before GWs are emitted, our estimates of the stochastic background only include BHs for which the expected number of coalescences during the instability time scale is $N_m < 1$, and for which $\tau_{\text{inst}} < \Delta t$ (which ensures that the instability time scale is shorter than the merger and accretion time scales).

The total energy emitted by the boson cloud during the signal duration $\Delta t$ can be estimated by integrating the GW energy flux given by Eq. (28). Using Eq. (29) we have

$$ \frac{dE_{\text{GW}}}{dt} = \frac{d\tilde{E} M_S^2}{dt M_f^2} = \frac{M_S^{\text{max}} \tau_{\text{GW}}}{(t + \tau_{\text{GW}})^2}, \quad (66) $$

and by integrating over a time $\Delta t$ we get

$$ E_{\text{GW}} = \int_{0}^{\Delta t} dt \frac{dE_{\text{GW}}}{dt} = \frac{M_{S}^{\text{max}} \Delta t}{\Delta t + \tau_{\text{GW}}}. \quad (67) $$

As shown in [79], the order of magnitude of the stochastic background can be estimated by computing the mass fraction of an isolated BH that is emitted by the boson cloud through GWs. This can be defined as

$$ f_{\text{ax}} = \frac{E_{\text{GW}}}{M_i}, \quad (68) $$

where we recall that $M_i$ is the initial mass of the BH. In Fig. 8 we show the average $f_{\text{ax}}$, weighted by the BH population, for our most optimistic models. In the LIGO and LISA band $f_{\text{ax}}$ can be order $\mathcal{O}(1\%)$, leading to a very large stochastic background [79].

Note that Eq. (64) cannot be applied to Galactic BHs which emit in the LIGO band, because it implicitly assumes that the number density of sources, $d^2\dot{n}/dMd\chi$, is homogeneous and isotropic. That assumption is clearly invalid for Galactic BHs [cf. Eq. (58)]. However, in this case we can simply sum the GW densities produced by the
<figure><img src="image_9.png" /><figcaption>FIG. 8. Average fraction of mass of an isolated BH emitted by the bosonic cloud for the optimistic models.</figcaption></figure>

Galactic BH population at the position of the detector. These densities are simply given by $\rho_{gw} = \dot{E}/(4\pi r^2) = (5/4)\pi f_s^2 h^2$ [cf. Eq. (42)], $r$ being the distance from the source to the detector. (Note that we neglect redshift and cosmological effects, since those are negligible inside the Galaxy.) Therefore, the GW energy density per (logarithmic) unit of frequency coming from each BH in the Galaxy is simply $d\rho_{gw}/d\ln f \approx (5/4)\pi f_s^2 h^2 \delta(\ln f - \ln f_s)$, where the Dirac delta is “spread out” over a frequency window of size $\sim \max[1/\Delta t, 1/T_{obs}]$ to account for the finite duration of the signal and the finite frequency resolution of the detector. Therefore, the contribution to the stochastic background from the population of Galactic BHs can be written as

$$ \Omega_{gw}(f) = \frac{1}{\rho_c} \int dM dV \frac{dn_{MW}}{dM} \frac{\Delta t}{t_0} \frac{d\rho_{gw}}{d\ln f}. \quad (69) $$

Here $dV$ denotes a volume integration over the Galaxy, and $\Delta t/t_0$ is again a duty cycle (i.e., we assume that Galactic BHs emit via boson condensates only once in their cosmic history).

To compute the SNR for the stochastic background we use

$$ \rho_{\text{stoch}} = \sqrt{T_{\text{obs}} \int_{f_{\text{min}}}^{f_{\text{max}}} df \frac{\Omega_{\text{GW}}^2}{\Omega_{\text{sens}}^2}}. \quad (70) $$

For LISA we have [121]

$$ \Omega_{\text{sens}} = S_h(f) \frac{2\pi^2}{3H_0^2} f^3, \quad (71) $$

while for LIGO [122]

$$ \Omega_{\text{sens}} = \frac{S_h(f)}{\sqrt{2\Gamma_{IJ}(f)}} \frac{2\pi^2}{3H_0^2} f^3, \quad (72) $$

where LIGO’s noise power spectral density $S_h(f)$ is assumed to be the same for both Livingston and Hanford, and $\Gamma_{IJ}$ is the overlap reduction function as defined in [123]. Notice the $1/\sqrt{2}$ factor in $\Omega_{sens}$ for LIGO compared to LISA, due to the use of data from two detectors instead of one.

As shown in Fig. 2 of [79], the SNR for this stochastic signal can be very high. Since the galactic background only contributes to the full spectrum in a very narrow frequency window around $f_s$, the contribution of the extragalactic background to the SNR largely dominates. When computing the background for LISA we assumed the semianalytic accretion model (C.1). Considering the most pessimistic accretion model (C.3) lowers the maximum SNR by at most a factor two.

## VI. EXCLUDING OR MEASURING BOSON MASSES THROUGH LISA BLACK HOLE SPIN MEASUREMENTS

So far we have focused on the direct detection of GWs from bosonic condensates. However it is also possible to infer the existence of light bosons in an indirect way. As shown in Fig. 1, the existence of a light boson would lead to the absence of BHs with spin above the corresponding superradiant instability window (i.e., there would be holes in the BH mass-spin “Regge plane” [17]). In this section we show that LISA measurements of the spins of merging massive BHs can be used to either rule out bosonic fields in the mass range $[4.5 \times 10^{-19}, 7.1 \times 10^{-13}]$ eV, or even more excitingly (if fields in the mass range $[10^{-17}, 10^{-13}]$ eV exist in nature) to measure their mass with percent accuracy.

In principle we could carry out a similar analysis using astrophysical models for stellar-mass BH binary mergers detectable by Advanced LIGO or third-generation Earth-based detectors. However, spin magnitude measurements for the components of a merging BH binaries are expected to be poor ($\Delta\chi \sim 0.3$ at best) even with third-generation detectors [40,41]. In addition, the mass range of BHs detectable by LIGO or future Earth-based interferometers overlaps in mass with existing spin estimates from low-mass x-ray binaries (see [38,124–126] for reviews of current BH spin estimates). In summary, we focus on LISA for two main reasons:

(i) LISA allows for percent-level determinations of massive BH spins (see e.g. Fig. 9 of [39]).

(ii) In comparison with current electromagnetic estimates of massive BH spins, which can be used to exclude boson masses in the range $[10^{-20}, 10^{-17}]$ eV (see e.g. [21,32]), LISA BH spin measurements can probe lower BH masses; therefore, depending on the details of massive BH formation models, they can exclude (or measure) boson masses all the way up to $m_s \sim 7 \times 10^{-13}$ eV.
One of our main tasks in this context is to determine whether LISA observations can distinguish between two models: one where a massive boson exists (depleting the corresponding instability region in the BH Regge plane) and a “standard” model where no depletion occurs. This is a standard Bayesian model selection problem (see e.g. [127–129] for previous applications of model selection to LISA observations of massive BH binaries).

We simulate massive BH binary catalogs corresponding to the three astrophysical models described in Sec. III B (popIII, Q3, Q3-nod) and seven values of $m_s$ in total, one for each decade in the boson mass range $m_s \in [10^{-19}, 10^{-13}]$ eV.

To simulate the loss of mass and angular momentum for each BH in the catalogs we compute the final angular momentum $J_f$ and mass $M_f$ according to Eqs. (25) and (26), with azimuthal number $1 \le m \le 4$ and frequency given by (8) with $l = m$ and $n = 0$. Approximating $\omega_R \approx \mu$ in Eq. (8) (which is strictly valid if $M\mu \ll 1$, but which is a good approximation even for $M\mu$ of order unity) we get

$$ \chi_f = \frac{4M_i\mu(m - M_i\chi_i)}{m^2}, \qquad (73) $$

$$ M_f = \frac{m - \sqrt{m^2 - 4mM_i\mu\chi_f + 4M_i^2\mu^2\chi_i\chi_f}}{2\mu\chi_f}. \qquad (74) $$

We migrate BHs in the Regge plane if the age of the Universe $t(z)$ at the merger redshift is larger than the instability time scale ($t(z) > \tau_{inst} = 1/\omega_I$) and if the spin is higher than a threshold $\chi_{max}(M, m_s)$ set by Eq. (25). This migration causes BHs in the catalog to accumulate along the critical line $\chi_{max}(M, m_s)$ in the Regge plane. An example of this accumulation can be seen in Fig. 9.

To compare two models $M_1$ and $M_2$ given a set of observations (i.e., a data set $D$), we can use Bayes’ theorem. The probability of model $M_i$ given the observations is

$$ P(M_i|D) = \frac{P(D|M_i)P(M_i)}{P(D)}, \qquad (75) $$

where $P(M_i)$ is the prior on model $M_i$, $P(D|M_i)$ is the likelihood of the data given the model, and $P(D)$ is an overall probability of observing the data $D$. Given a likelihood function for each model, we can then compute the odds ratio between the two models:

$$ \mathcal{O}(M_1/M_2) = \frac{P(M_1|D)}{P(M_2|D)} = \frac{P(D|M_1)P(M_1)}{P(D|M_2)P(M_2)}. \qquad (76) $$

A value of the odds ratio larger than one favors model $M_1$, while a value of the odds ratio lower than one favors model $M_2$. When $P(M_1) = P(M_2)$ the last factor on the right-hand side simplifies, and the odds ratio is just the ratio of

<figure><img src="image_10.png" /><figcaption>FIG. 9. Example of a two-year simulation of massive BHs as observed by LISA assuming the Q3-nod model in the presence of a boson of mass $m_s = 10^{-16}$ eV. Each blue circle corresponds to the mass and spin of one component of an observed BH binary. The brown line corresponds to the maximum allowed spin $\chi_{max}(M, m_s)$ for the given boson mass. This curve is shaped like a sawtooth because different $m$-harmonics are more important for different BH masses. In this particular instance, LISA measurements from the simulated data would lead to a measured boson mass $0.88 \times 10^{-16}$ eV $< m_s^m < 1.35 \times 10^{-16}$ eV.</figcaption></figure>

the likelihood of the data in both models (also known as the “Bayes factor”).

We construct a likelihood function for BHs in the Regge plane for two models: one with no ultralight boson, and one with an ultralight boson of mass $m_s$. To avoid a possible bias toward high spins in the astrophysical models (see e.g. [76]) we choose the simplest likelihood function in the absence of bosons: $\mathcal{L}_0(M, \chi) = 1$. In the presence of bosons, we set the likelihood $\mathcal{L}_{m_s}(M, \chi)$ to unity if $\chi \le \chi_{max}(M, m_s)$, and we set it to zero otherwise. We add to this likelihood a Gaussian centered on the threshold $\chi_{max}(M, m_s)$ with width $\sigma_\chi=0.05$, with a prefactor $1 - \chi_{max}$ in front of it. This factor represents the fraction of BHs with spins higher than the threshold that have migrated out of the exclusion region to accumulate on the threshold line, under the simplifying assumption that they migrate in the $\chi$ direction only (i.e., we neglect the relatively small variations in the BH mass). In summary, the likelihood $\mathcal{L}_{m_s}(M, \chi)$ in the presence of a boson of mass $m_s$ is defined by

$$ \mathcal{L}_{m_s}(M, \chi) = \begin{cases} 1, & \chi_{\max}(M, m_s) = 1 \\ 1 + G(\chi, 0.05), & \chi < \chi_{\max}(M, m_s) < 1, \\ G(\chi, 0.05), & \chi_{\max}(M, m_s) < \chi < 1 \end{cases} \qquad (77) $$

$$ G(\chi, \sigma) = \frac{1 - \chi_{\max}}{\sqrt{2\pi}\sigma} \exp \left[ - \frac{(\chi - \chi_{\max})^2}{2\sigma^2} \right]. \qquad (78) $$
The prefactor in front of the Gaussian ensures that the two likelihoods $\mathcal{L}_0(M, \chi)$ and $\mathcal{L}_{m_s}(M, \chi)$ have the same “weight,” in the sense that the integral $\int \mathcal{L}_M dM d\chi$ is independent of the model (so the presence or absence of an ultralight boson have, a priori, the same probability). This choice for the likelihood functions assures that the computation of the odds ratio is agnostic about the underlying astrophysical model.

As stated earlier, the spin threshold $\chi_{\text{max}}(M, m_s)$ is given by Eq. (25). In practice this criterion is slightly complicated by the fact that the range of affected BH masses depends on the time available for each system to radiate, which in turn depends on the redshift. For simplicity we compute the spin limit using a constant instability time scale of 500 Myrs (approximately the age of the Universe at redshift $z=10$), setting $\omega_R = \mu$ in Eq. (25). The choice of this time scale is conservative in the sense that the exclusion region is smaller than it would have been if we had chosen longer time scales. Indeed, our choice reduces the likelihood discrepancy for low redshift BHs that will have migrated to the threshold line, but would not have had the time to do so had they merged at higher redshifts. For illustration, Fig. 9 shows the distribution of BH masses and spins for one realization of a two-year catalog with $m_s = 10^{-16}$ eV, along with the corresponding spin threshold $\chi_{\text{max}}(M, m_s)$.

We simulate LISA observations of these catalogs using a Fisher-matrix analysis similar to the study presented in [39], using the updated LISA noise PSD of [11]. In addition to instrumental noise, we also include the boson mass-dependent confusion noise coming from superradiant BH instabilities shown in Fig. 1 of [79]. For each detectable binary (where detectability is defined as $\rho > 10)^6$ we approximate the recovered distribution for each binary BH component by a bivariate Gaussian centered on the true values $(\bar{M}_i, \bar{\chi}_i)$, with spread given by the two-dimensional inverse of the covariance matrix $\Gamma = \Sigma^{-1}$:

$$p_{\text{obs}}(M_i, \chi_i) = \frac{\sqrt{|\Gamma|}}{2\pi} \exp\left\{ -\frac{1}{2} \left[ \begin{aligned} & \Gamma_{M_i M_i} (M_i - \bar{M}_i)^2 \\ & + \Gamma_{\chi_i \chi_i} (\chi_i - \bar{\chi}_i)^2 \\ & + 2\Gamma_{M_i \chi_i} (M_i - \bar{M}_i)(\chi_i - \bar{\chi}_i) \end{aligned} \right] \right\}. \quad (79)$$

One problem is that GW observations can measure the reshifted mass $M_z = (1+z)M$, rather than the BH mass in the source frame $M$. Lensing effects will induce an extra uncertainty on the distance to the source of typical size $\sigma_{D_L}^{\text{lens}}(z)$, and through the redshift-distance relation $D_L(z)$ an extra uncertainty on the redshift of size $\sigma_z^{\text{lens}}(z)$. We include

the effects of lensing by adjusting the observed distribution $p_{\text{obs}}(M_i, \chi_i)$ along the mass direction. We estimate the typical extra error on the mass due to lensing as

$$\frac{\sigma_M^{\text{lens}}(z)}{M} = \frac{\sigma_z^{\text{lens}}(z)}{1+z} = \frac{dz}{dD_L}(z) \frac{\sigma_{D_L}^{\text{lens}}(z)}{1+z}. \quad (80)$$

where the luminosity distance error as a function of redshift can be estimated by the approximate relation [44,130]

$$\sigma_{D_L}^{\text{lens}}(z) = D_L(z) \times 0.066\{4[1 - (1+z)^{-1/4}]\}^{9/5}. \quad (81)$$

At this stage we can compute the likelihood of an observed BH for each model $M$ by integrating the product

$$\mathcal{L}(i|M) = \int p_{\text{obs}}(M_i, \chi_i) \mathcal{L}_M(M_i, \chi_i) dM_i d\chi_i, \quad (82)$$

where the index $i$ labels the observed BH. In the absence of ultralight bosons we get $\mathcal{L}(i|M_0) = 1$, and in the presence of bosons we use Monte Carlo methods to compute $\mathcal{L}(i|m_s)$. In practice we generate a set of random points in the Regge plane $(M_k, \chi_k)$ distributed according to $p_{\text{obs}}(M_i, \chi_i)$, with an extra (spin-independent) jump in the mass direction due to lensing, which we assume to be Gaussian distributed with zero mean and standard deviation $\sigma_M^{\text{lens}}(z)$. The integral is then approximated by

$$\mathcal{L}(i|m_s) \approx \frac{1}{N} \prod_{k=1}^{N} \mathcal{L}_{m_s}(M_k, \chi_k). \quad (83)$$

The integration with respect to mass and spin in Eq. (82) tends to suppress the effect on the odds ratio of potential observations in the exclusion region that would favor high spins with low confidence. As one can see from Eq. (79), if the measurement error on the spin is significant, Eq. (82) will show a significant overlap between the two factors inside the integral, even if the observed spin is higher than the threshold.

Using this method we can simulate a set of LISA observations $D$ and compute its likelihood for model $M$ as

$$\mathcal{L}(D|M) = \prod_i \mathcal{L}(i|M), \quad (84)$$

where the product is taken over all components of a binary observed with SNR $\rho > 10$. Then, assuming no prior preference, we compute the odds ratio between a model with boson mass $m_s$ and a model without bosons:

$$\mathcal{O}(m_s/M_0) = \mathcal{L}(D|m_s). \quad (85)$$

We simulated observations in the absence of ultralight bosons and in the presence of an ultralight boson with

<footnote>$^6$Note that this threshold is slightly different from that used elsewhere in the paper ($\rho = 8$, though that was for boson-condensate sources). Still, the results hardly depend on this choice, since barely detectable events ($\rho \sim 8-10$) have anyhow very poor spin determinations.</footnote>
<table><caption>TABLE III. Median minimum and maximum boson mass excluded by LISA for different BH evolution models (popIII, Q3, Q3-nod) and observation times $T_{\text{obs}}$.</caption><thead><tr><th>Model</th><th><i>T</i><sub>obs</sub> [yr]</th><th>Min [eV]</th><th>Max [eV]</th></tr></thead><tbody><tr><td rowspan="4">popIII</td><td>0.5</td><td>4.7 × 10<sup>-17</sup></td><td>4.7 × 10<sup>-14</sup></td></tr><tr><td>1</td><td>8.2 × 10<sup>-18</sup></td><td>8.9 × 10<sup>-14</sup></td></tr><tr><td>2</td><td>6.9 × 10<sup>-18</sup></td><td>1.3 × 10<sup>-13</sup></td></tr><tr><td>4</td><td>4.5 × 10<sup>-18</sup></td><td>1.6 × 10<sup>-13</sup></td></tr><tr><td rowspan="4">Q3</td><td>0.5</td><td>—</td><td>—</td></tr><tr><td>1</td><td>9.4 × 10<sup>-18</sup></td><td>1.9 × 10<sup>-15</sup></td></tr><tr><td>2</td><td>6.9 × 10<sup>-18</sup></td><td>7.5 × 10<sup>-15</sup></td></tr><tr><td>4</td><td>4.1 × 10<sup>-18</sup></td><td>8 × 10<sup>-15</sup></td></tr><tr><td rowspan="4">Q3-nod</td><td>0.5</td><td>6.9 × 10<sup>-18</sup></td><td>3.6 × 10<sup>-15</sup></td></tr><tr><td>1</td><td>4.5 × 10<sup>-18</sup></td><td>6.7 × 10<sup>-15</sup></td></tr><tr><td>2</td><td>1.8 × 10<sup>-18</sup></td><td>1 × 10<sup>-14</sup></td></tr><tr><td>4</td><td>1 × 10<sup>-18</sup></td><td>2.3 × 10<sup>-14</sup></td></tr></tbody></table>

seven possible values of $m_s$ (one for each decade in the boson mass range $m_s \in [10^{-19}, 10^{-13}]$ eV). For each boson mass and for the model without bosons, we simulated a set of 21 realizations of the LISA mission considering the three astrophysical models (popIII, Q3, and Q3-nod) and four choices for the observation time (6 months, 1 year, 2 years and 4 years), corresponding to a total of 252 simulations per model.

In the absence of an ultralight boson, we identify as excluded the range of masses where the odds ratio $\mathcal{O}$ satisfies $\log[\mathcal{O}(m_s/M_0)] < -4.5$. This criterion corresponds to rejecting the presence of the given boson mass at 3-$\sigma$ confidence level. This requirement to exclude a boson of a given mass corresponds to a false alarm rate of ~10% for a four-year mission in the popIII model, and less than 5% in the other models: the maximum odds ratio incorrectly favoring the presence of an ultralight boson in the 84 realizations where we assumed its absence was $\log(O_{\text{max}}) = 5.2$ for the popIII model, 1.1 for the Q3 model, and 2.8 for the Q3-nod model. In the popIII case, a maximum odds ratio of 4.5 was exceeded twice. The median range of boson masses excluded in our simulations is summarized in Table III and Fig. 10. As expected, in our light-seed (popIII) model the excluded boson masses are higher than in the heavy-seed models Q3 and Q3-nod, because the observed BH masses are generally lower in light-seed scenarios. The Q3-nod model allows us to set more stringent bounds than the Q3 model, because the merger rate is higher when there are no delays between galaxy mergers and BH mergers. Furthermore, the Q3 model failed to allow for a boson mass exclusion after six months of observations in 12 of the 21 simulations due to its low merger rate. For any astrophysical models among the three we considered, four years of LISA observations would allow us to exclude boson masses in the range $4.1 \times 10^{-18}$ eV to $8 \times 10^{-15}$ eV.

<figure><img src="image_11.png" /><figcaption>FIG. 10. Median minimum and maximum boson mass excluded by LISA for different observation times $T_{\text{obs}}$ and BH evolution models (red, solid line: popIII; green, dotted line: Q3; blue, dashed line: Q3-nod). Due to the reduced merger rate in the Q3 model, limits on the boson mass could be put in more than half of the simulations only after one year of observations.</figcaption></figure>

It is also interesting to address the following question: in the presence of an ultralight boson, could it be detected? And if so, what is the accuracy with which we could determine its mass? To answer the first question we identify the mass range where $\log[\mathcal{O}(m_s)/O_{\text{max}}] \ge -4.5$, again corresponding to a 3-$\sigma$ confidence level, and then use the simulated events to determine the accuracy with which $m_s$ can be determined. Our results are summarized in Table IV. We do not show results for masses where four years of observations were not enough to claim a boson detection. In marginal detections ($\log(O_{\text{max}}) \lesssim 10$), only the order of magnitude of the boson mass could be inferred.

For the light-seed popIII model, boson masses in the range $[10^{-16}, 10^{-13}]$ eV could be confidently detected after four years of observations with measurement errors in $m_s$ of 5–10%. Model Q3-nod allows for the confident detection of a boson in the mass range $[10^{-17}, 10^{-14}]$ eV with mass measurement errors of 5–15%, while the less optimistic model Q3 only allows detections for bosons with mass in the range $[10^{-16}, 10^{-15}]$ eV, with mass measurement errors of ~40%.

We remark that the biases in the recovered boson masses are sometimes comparable to the corresponding measurement accuracies: in low-mass (high-mass) seed models we tend to overestimate (underestimate) the boson mass. It is likely that this bias could be reduced with better modeling of the relevant physics—e.g. by evolving Eqs. (24) numerically for each BH from formation until merger—or with a more careful choice of the likelihood function, e.g. by taking the observed redshift of the system into account in the definition of the threshold line in Eq. (77), i.e., in the likelihood in the presence of bosons. A more detailed analysis of systematic and statistical errors in recovering the boson masses is an interesting topic for future work.
<table><caption>TABLE IV. Median measured boson mass $m_s^m = \kappa m_s$ and median maximum log likelihood $L = \log O_{\text{max}}$ for different BH evolution models, observation times $T_{\text{obs}}$, and “true” boson masses $m_s$.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">$T_{\text{obs}}$ [yr]<br>$m_s$ [eV]</th><th colspan="2">0.5</th><th colspan="2">1</th><th colspan="2">2</th><th colspan="2">4</th></tr><tr><th>$\kappa$</th><th>$L$</th><th>$\kappa$</th><th>$L$</th><th>$\kappa$</th><th>$L$</th><th>$\kappa$</th><th>$L$</th></tr></thead><tbody><tr><td rowspan="4">popIII</td><td>$10^{-16}$</td><td>—</td><td>2.9</td><td>—</td><td>4.1</td><td>1.06 ± 0.25</td><td>13</td><td>1.07 ± 0.12</td><td>28</td></tr><tr><td>$10^{-15}$</td><td>1 ± 0.4</td><td>7.9</td><td>1.05 ± 0.21</td><td>14</td><td>1.06 ± 0.11</td><td>39</td><td>1.08 ± 0.06</td><td>90</td></tr><tr><td>$10^{-14}$</td><td>1 ± 0.6</td><td>5.4</td><td>1.02 ± 0.15</td><td>12</td><td>1.05 ± 0.1</td><td>31</td><td>1.06 ± 0.06</td><td>81</td></tr><tr><td>$10^{-13}$</td><td>—</td><td>0.64</td><td>—</td><td>1.7</td><td>1 ± 0.15</td><td>8.6</td><td>1.02 ± 0.1</td><td>26</td></tr><tr><td rowspan="2">Q3</td><td>$10^{-16}$</td><td>—</td><td>0.91</td><td>—</td><td>3.2</td><td>—</td><td>4.5</td><td>1 ± 0.4</td><td>9.7</td></tr><tr><td>$10^{-15}$</td><td>—</td><td>0</td><td>—</td><td>1.9</td><td>—</td><td>3.6</td><td>1 ± 0.4</td><td>6.8</td></tr><tr><td rowspan="4">Q3-nod</td><td>$10^{-17}$</td><td>—</td><td>2.9</td><td>1 ± 0.23</td><td>6.5</td><td>1.03 ± 0.19</td><td>13</td><td>1.02 ± 0.13</td><td>25</td></tr><tr><td>$10^{-16}$</td><td>1 ± 0.4</td><td>17</td><td>0.99 ± 0.15</td><td>47</td><td>1 ± 0.08</td><td>98</td><td>0.97 ± 0.06</td><td>200</td></tr><tr><td>$10^{-15}$</td><td>1 ± 0.5</td><td>11</td><td>0.94 ± 0.18</td><td>28</td><td>0.95 ± 0.1</td><td>65</td><td>0.98 ± 0.07</td><td>140</td></tr><tr><td>$10^{-14}$</td><td>—</td><td>1.6</td><td>—</td><td>4.2</td><td>0.98 ± 0.21</td><td>14</td><td>0.98 ± 0.13</td><td>27</td></tr></tbody></table>

## VII. CONCLUSIONS AND OUTLOOK

In this work and in the companion paper [79] we assess the detectability of light-boson condensates around BHs with GW interferometers combining the best available estimates for GW emission from these systems, state-of-the-art astrophysical BH population models, and relatively realistic GW data analysis techniques.

For both Advanced LIGO and LISA, we find that the most stringent constraints on the boson mass $m_s$ should come from the stochastic background produced by the superposition of unresolved GW signals from BH-boson condensate systems. We show that this background should be detectable by Advanced LIGO for $m_s \in [2 \times 10^{-13}, 10^{-12}]$ eV, and by LISA for $m_s \in [5 \times 10^{-19}, 5 \times 10^{-16}]$ eV. We also find that existing constraints on the stochastic background from Advanced LIGO’s O1 run may already rule out a range of boson masses in the Advanced LIGO window. Although the precise constrained regions depend on the astrophysical model, the order of magnitude of the stochastic background is robust with respect to astrophysical uncertainties, as shown in Fig. 2 of [79].

Our results indicate that ~15–200 resolvable sources should be detectable by Advanced LIGO for scalar field masses $m_s \sim 3 \times 10^{-13}$ eV, while LISA should be able to resolve ~5–40 sources for $m_s \sim 10^{-17}$ eV. Moreover, LISA measurements of BH spins may either determine $m_s \in [10^{-17}, 10^{-13}]$ eV to within 10% accuracy, or rule out boson masses in the range $m_s \in [10^{-18}, 1.6 \times 10^{-13}]$ eV.

We anticipate that pulsar-timing arrays [131–135], though sensitive to the stochastic GW background in the nHz band, may not set stringent constraints on the masses of ultralight bosons. The reason lies in the very large instability and gravitational radiation time scales for bosons masses in the nHz band and in the paucity of massive BHs with masses $M \gtrsim 10^{10} M_{\odot}$ [136,137], which would be required to produce a significant background

from BH-boson condensates. Conversely, an interferometer like DECIGO [138] would allow one to put constraints on boson masses $m_s \sim 10^{-14}$ eV.

Some of our conclusions differ from previous work on this topic by Arvanitaki et al. [19,20], which neglected the stochastic background from boson condensates in the LISA and in the LIGO band, focusing on resolved events. This had the two-fold effect of (i) missing the strong constraints (summarized above) from existing and projected stochastic background limits, and (ii) missing the “self-confusion” problem, i.e., the fact that the stochastic background itself is a confusion noise (similar to the familiar white dwarf confusion noise in the LISA band), impairing the detectability of individual sources.

Another important difference with respect to Arvanitaki et al. [19,20] lies in our astrophysical models. Refs. [19,20] focused on Galactic BHs as resolvable LIGO sources. This is probably the main reason why they overlooked the presence of a significant stochastic background, which is mostly produced by extragalactic BHs. Likewise, the lower LISA event rates found by [19] (in spite of their neglecting the aforementioned confusion noise from the background) seem to be due to their simplified (and overly pessimistic) models for the massive BH population. For example, Ref. [19] considered the chaotic accretion model of [139], where BHs with large spins are unlikely. Such models are either disfavored or ruled out (depending on the assumed spin distribution) by iron $K\alpha$ line data [76].

Finally, our analysis of the statistical error affecting GW measurements of BH spins in the LISA band and our use of Bayesian model selection techniques (while far from realistic) are a step forward with respect to the estimates of [20], and they lead to one of the most remarkable conclusions of our work. As shown schematically in Fig. 1, LISA could either rule out light bosons in the mass range $[4 \times 10^{-18}, 10^{-14}]$ eV, or measure $m_s$ with ten percent
accuracy if particles in the mass range [10⁻¹⁷, 10⁻¹³] eV exist in Nature.

## ACKNOWLEDGMENTS

We thank Asimina Arvanitaki, Masha Baryakhtar and Robert Lasenby for useful comments. S. Ghosh and E. Berti are supported by NSF Grants No. PHY-1607130, AST-1716715 and by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Programme. V. Cardoso acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant "Matter and strong-field gravity: New frontiers in Einstein's theory" grant agreement no. MaGRaTh-646597. The work of I. Dvorkin has been done within the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d'avenir under the reference ANR-11-IDEX-0004-02. Research at Perimeter

Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690904 and from FCT-Portugal through the project IF/00293/2013. The authors would like to acknowledge networking support by the COST Action CA16104. This work has made use of the Horizon Cluster, hosted by the Institut d'Astrophysique de Paris. We thank Stephane Rouberol for running smoothly this cluster for us. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by Sérgio Almeida at CENTRA/IST. Computations were performed at the cluster "Baltasar-Sete-Sóis", and supported by the MaGRaTh-646597 ERC Consolidator Grant. This work was supported by the French Centre National d'Études Spatiales (CNES).

[1] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. **116**, 061102 (2016).

[2] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. **116**, 241103 (2016).

[3] B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev. Lett. **118**, 221101 (2017).

[4] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. X **6**, 041015 (2016).

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

[6] V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, Phys. Rev. D **95**, 084014 (2017).

[7] V. Cardoso, E. Franzin, and P. Pani, Phys. Rev. Lett. **116**, 171101 (2016); **117**, 089902(E) (2016).

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

[9] C. Chirenti and L. Rezzolla, Phys. Rev. D **94**, 084016 (2016).

[10] P. Amaro-Seoane et al., GW Notes **6**, 4 (2013).

[11] H. Audley et al., arXiv:1702.00786.

[12] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D **81**, 123530 (2010).

[13] R. Essig et al., arXiv:1311.0029.

[14] D. J. E. Marsh, Phys. Rep. **643**, 1 (2016).

[15] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D **95**, 043541 (2017).

[16] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. **906**, pp. 1 (2015).

[17] A. Arvanitaki and S. Dubovsky, Phys. Rev. D **83**, 044026 (2011).

[18] R. Brito, V. Cardoso, and P. Pani, Classical Quantum Gravity **32**, 134001 (2015).

[19] A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D **91**, 084011 (2015).

[20] A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky, and R. Lasenby, Phys. Rev. D **95**, 043001 (2017).

[21] M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D **96**, 035019 (2017).

[22] P. S. B. Dev, M. Lindner, and S. Ohmer, arXiv:1609.03939.

[23] R. Penrose, Riv. Nuovo Cimento **1**, 252 (1969); Gen. Relativ. Gravit. **34**, 1141 (2002).

[24] Y. B. Zel'dovich, Pis'ma Zh. Eksp. Teor. Fiz. **14**, 270 (1971) [JETP Lett. **14**, 180 (1971)].

[25] Y. B. Zel'dovich, Zh. Eksp. Teor. Fiz. **62**, 2076 (1972) [Sov. Phys. JETP **35**, 1085 (1972)].

[26] W. H. Press and S. A. Teukolsky, Nature (London) **238**, 211 (1972).

[27] V. Cardoso, O. J. C. Dias, J. P. S. Lemos, and S. Yoshida, Phys. Rev. D **70**, 044039 (2004).

[28] S. L. Detweiler, Phys. Rev. D **22**, 2323 (1980).

[29] T. J. M. Zouros and D. M. Eardley, Ann. Phys. (N.Y.) **118**, 139 (1979).

[30] V. Cardoso and S. Yoshida, J. High Energy Phys. 07 (2005) 009.

[31] S. R. Dolan, Phys. Rev. D **76**, 084001 (2007).

[32] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. Lett. **109**, 131102 (2012).

[33] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. D **86**, 104017 (2012).

[34] W. E. East and F. Pretorius, Phys. Rev. Lett. **119**, 041101 (2017).

[35] W. E. East, Phys. Rev. D **96**, 024004 (2017).

[36] R. Brito, V. Cardoso, and P. Pani, Phys. Rev. D **88**, 023514 (2013).
[37] L. Brenneman, C. Reynolds, M. Nowak, R. Reis, M. Trippe et al., Astrophys. J. **736**, 103 (2011).

[38] M. Middleton, Black Hole Spin: Theory and Observation, in *Astrophysics of Black Holes*, edited by C. Bambi, Astrophysics and Space Science Library Vol. 440 (Springer, Berlin, Heidelberg, 2016), pp. 99–151.

[39] A. Klein et al., Phys. Rev. D **93**, 024003 (2016).

[40] S. Vitale, R. Lynch, V. Raymond, R. Sturani, J. Veitch, and P. Graff, Phys. Rev. D **95**, 064053 (2017).

[41] S. Vitale and M. Evans, Phys. Rev. D **95**, 064052 (2017).

[42] H. Yoshino and H. Kodama, Prog. Theor. Exp. Phys. **2014**, 043E02 (2014).

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

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

[45] E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. **117**, 101102 (2016).

[46] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D **95**, 103012 (2017).

[47] J. R. Gair, S. Babak, A. Sesana, P. Amaro-Seoane, E. Barausse, C. P. L. Berry, E. Berti, and C. Sopuerta, J. Phys. Conf. Ser. **840**, 012021 (2017).

[48] H. Okawa, H. Witek, and V. Cardoso, Phys. Rev. D **89**, 104032 (2014).

[49] M. Zilhão, H. Witek, and V. Cardoso, Classical Quantum Gravity **32**, 234003 (2015).

[50] N. Sanchis-Gual, J. C. Degollado, P. J. Montero, J. A. Font, and C. Herdeiro, Phys. Rev. Lett. **116**, 141101 (2016).

[51] P. Bosch, S. R. Green, and L. Lehner, Phys. Rev. Lett. **116**, 141102 (2016).

[52] D. Christodoulou and R. Ruffini, Phys. Rev. D **4**, 3552 (1971).

[53] C. A. R. Herdeiro and E. Radu, arXiv:1706.06597.

[54] E. Berti, V. Cardoso, and M. Casals, Phys. Rev. D **73**, 024013 (2006).

[55] E. Berti, V. Cardoso, and A. O. Starinets, Classical Quantum Gravity **26**, 163001 (2009).

[56] S. A. Teukolsky, Astrophys. J. **185**, 635 (1973).

[57] M. Sasaki and H. Tagoshi, Living Rev. Relativ. **6**, 6 (2003).

[58] E. Poisson and M. Sasaki, Phys. Rev. D **51**, 5753 (1995).

[59] S. A. Teukolsky and W. H. Press, Astrophys. J. **193**, 443 (1974).

[60] N. Sanchis-Gual, J. C. Degollado, C. Herdeiro, J. A. Font, and P. J. Montero, Phys. Rev. D **94**, 044061 (2016).

[61] N. Sanchis-Gual, J. C. Degollado, J. A. Font, C. Herdeiro, and E. Radu, Classical Quantum Gravity **34**, 165001 (2017).

[62] C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. **112**, 221101 (2014).

[63] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D **24**, 1542014 (2015).

[64] C. Herdeiro, E. Radu, and H. Runarsson, Classical Quantum Gravity **33**, 154001 (2016).

[65] B. P. Abbott et al. (Virgo, LIGO Scientific), arXiv: 1707.02667.

[66] C. Palomba (VIRGO, LIGO Scientific), arXiv:1201.3176.

[67] A. R. King and U. Kolb, Mon. Not. R. Astron. Soc. **305**, 654 (1999).

[68] M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel, and R. O'Shaughnessy, Astrophys. J. **759**, 52 (2012).

[69] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O'Shaughnessy, Astrophys. J. **779**, 72 (2013).

[70] M. Dominik, E. Berti, R. O'Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Astrophys. J. **806**, 263 (2015).

[71] K. Belczynski, S. Repetto, D. E. Holz, R. O'Shaughnessy, T. Bulik, E. Berti, C. Fryer, and M. Dominik, Astrophys. J. **819**, 108 (2016).

[72] K. Belczynski, D. E. Holz, T. Bulik, and R. O'Shaughnessy, Nature (London) **534**, 512 (2016).

[73] D. Gerosa and E. Berti, Phys. Rev. D **95**, 124046 (2017).

[74] M. Fishbach, D. E. Holz, and B. Farr, Astrophys. J. **840**, L24 (2017).

[75] E. Berti and M. Volonteri, Astrophys. J. **684**, 822 (2008).

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

[77] J. X. Wang and D. Merritt, Astrophys. J. **600**, 149 (2004).

[78] B. P. Abbott et al. (VIRGO, LIGO Scientific), Living Rev. Relativ. **19**, 1 (2016).

[79] R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, arXiv:1706.05097.

[80] A. Stroeer and A. Vecchio, Classical Quantum Gravity **23**, S809 (2006).

[81] A. J. Ruiter, K. Belczynski, M. Benacquista, S. L. Larson, and G. Williams, Astrophys. J. **717**, 1006 (2010).

[82] B. Behnke, M. A. Papa, and R. Prix, Phys. Rev. D **91**, 064007 (2015).

[83] B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. D **94**, 102002 (2016).

[84] C. Cutler, Phys. Rev. D **57**, 7089 (1998).

[85] E. Berti, A. Buonanno, and C. M. Will, Phys. Rev. D **71**, 084025 (2005).

[86] M. Maggiore, *Gravitational Waves. Vol. 1: Theory and Experiments*, Oxford Master Series in Physics (Oxford University Press, New York, 2007).

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

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

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

[90] 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).

[91] E. Barausse, F. Shankar, M. Bernardi, Y. Dubois, and R. K. Sheth, Mon. Not. R. Astron. Soc. **468**, 4782 (2017).

[92] I. Dvorkin and E. Barausse, Mon. Not. R. Astron. Soc. **470**, 4547 (2017).

[93] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, Astrophys. J. **178**, 347 (1972).

[94] F. Shankar, D. H. Weinberg, and J. Miralda-Escude', Mon. Not. R. Astron. Soc. **428**, 421 (2013).

[95] K. Pardo, A. D. Goulding, J. E. Greene, R. S. Somerville, E. Gallo, R. C. Hickox, B. P. Miller, A. E. Reines, and J. D. Silverman, Astrophys. J. **831**, 203 (2016).
[96] S. E. de Mink, M. Cantiello, N. Langer, O. R. Pols, I. Brott, and S. C. Yoon, Astron. Astrophys. **497**, 243 (2009).

[97] K. Belczynski, A. Buonanno, M. Cantiello, C. L. Fryer, D. E. Holz, I. Mandel, M. C. Miller, and M. Walczak, Astrophys. J. **789**, 120 (2014).

[98] I. Mandel and S. E. de Mink, Mon. Not. R. Astron. Soc. **458**, 2634 (2016).

[99] P. Marchant, N. Langer, P. Podsiadlowski, T. M. Tauris, and T. J. Moriya, Astron. Astrophys. **588**, A50 (2016).

[100] J. S. Vink, New Astron. **52**, 419 (2008).

[101] K. Belczynski, M. Dominik, T. Bulik, R. O'Shaughnessy, C. Fryer, and D. E. Holz, Astrophys. J. **715**, L138 (2010).

[102] F. Antonini and F. A. Rasio, Astrophys. J. **831**, 187 (2016).

[103] S. Chatterjee, C. L. Rodriguez, and F. A. Rasio, Astrophys. J. **834**, 68 (2017).

[104] C. L. Fryer, K. Belczynski, G. Wiktorowicz, M. Dominik, V. Kalogera, and D. E. Holz, Astrophys. J. **749**, 91 (2012).

[105] I. Dvorkin, E. Vangioni, J. Silk, J. P. Uzan, and K. A. Olive, Mon. Not. R. Astron. Soc. **461**, 3877 (2016).

[106] S. E. Woosley and T. A. Weaver, Astrophys. J. Suppl. Ser. **101**, 181 (1995).

[107] E. Vangioni, K. A. Olive, T. Prestegard, J. Silk, P. Petitjean, and V. Mandic, Mon. Not. R. Astron. Soc. **447**, 2575 (2015).

[108] R. J. Bouwens, G. D. Illingworth, P. A. Oesch, I. Labbé, M. Trenti, P. van Dokkum, M. Franx, M. Stiavelli, C. M. Carollo, D. Magee, and V. Gonzalez, Astrophys. J. **737**, 90 (2011).

[109] P. S. Behroozi, R. H. Wechsler, and C. Conroy, Astrophys. J. **770**, 57 (2013).

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

[111] D. Schaerer, Astron. Astrophys. **382**, 28 (2002).

[112] P. S. Behroozi, R. H. Wechsler, and C. Conroy, Astrophys. J. **762**, L31 (2013).

[113] X. Ma, P. F. Hopkins, C. A. Faucher-Giguère, N. Zolman, A. L. Muratov, D. Kereš, and E. Quataert, Mon. Not. R. Astron. Soc. **456**, 2140 (2016).

[114] L. Hernquist, Astrophys. J. **356**, 359 (1990).

[115] S. Shen, H. J. Mo, S. D. M. White, M. R. Blanton, G. Kauffmann, W. Voges, J. Brinkmann, and I. Csabai, Mon. Not. R. Astron. Soc. **343**, 978 (2003).

[116] C. Porcel, F. Garzon, J. Jimenez-Vicente, and E. Battaner, Astron. Astrophys. **330**, 136 (1998).

[117] T. Hartwig, M. Volonteri, V. Bromm, R. S. Klessen, E. Barausse, M. Magg, and A. Stacy, Mon. Not. R. Astron. Soc. **460**, L74 (2016).

[118] J. M. Bardeen, Nature (London) **226**, 64 (1970).

[119] J. Abadie et al. (VIRGO, LIGO Scientific), Classical Quantum Gravity **27**, 173001 (2010).

[120] E. S. Phinney, arXiv:astro-ph/0108028.

[121] N. J. Cornish and S. L. Larson, Classical Quantum Gravity **18**, 3473 (2001).

[122] B. Allen and J. D. Romano, Phys. Rev. D **59**, 102001 (1999).

[123] E. Thrane and J. D. Romano, Phys. Rev. D **88**, 124032 (2013).

[124] R. A. Remillard and J. E. McClintock, Annu. Rev. Astron. Astrophys. **44**, 49 (2006).

[125] C. S. Reynolds, Space Sci. Rev. **183**, 277 (2014).

[126] M. C. Miller and J. M. Miller, Phys. Rep. **548**, 1 (2015).

[127] J. R. Gair, A. Sesana, E. Berti, and M. Volonteri, Classical Quantum Gravity **28**, 094018 (2011).

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

[129] E. Berti, J. Gair, and A. Sesana, Phys. Rev. D **84**, 101501 (2011).

[130] C. M. Hirata, D. E. Holz, and C. Cutler, Phys. Rev. D **81**, 124046 (2010).

[131] M. Kramer and D. J. Champion, Classical Quantum Gravity **30**, 224009 (2013).

[132] G. Hobbs, Classical Quantum Gravity **30**, 224007 (2013).

[133] M. A. McLaughlin, Classical Quantum Gravity **30**, 224008 (2013).

[134] G. Hobbs et al., Classical Quantum Gravity **27**, 084013 (2010).

[135] R. N. Manchester and (IPTA), Classical Quantum Gravity **30**, 224010 (2013).

[136] P. Natarajan and E. Treister, Mon. Not. R. Astron. Soc. **393**, 838 (2009).

[137] F. Pacucci, P. Natarajan, and A. Ferrara, Astrophys. J. **835**, L36 (2017).

[138] S. Kawamura et al., Classical Quantum Gravity **23**, S125 (2006).

[139] A. R. King, J. E. Pringle, and J. A. Hofmann, Mon. Not. R. Astron. Soc. **385**, 1621 (2008).