Stochastic and resolvable gravitational waves from ultralight bosons

TL;DR

This paper calculates the stochastic gravitational-wave background from ultralight bosons around black holes, showing it could be detectable by LIGO and LISA, and may help constrain axion-like particles.

Contribution

First computation of the stochastic gravitational-wave background from superradiant bosonic condensates around black holes using numerical solutions and astrophysical models.

Findings

Background potentially observable by Advanced LIGO and LISA.

LIGO O1 data may already constrain certain axion masses.

Semicoherent searches could detect dozens of resolvable sources.

Abstract

Ultralight scalar fields around spinning black holes can trigger superradiant instabilities, forming a long-lived bosonic condensate outside the horizon. We use numerical solutions of the perturbed field equations and astrophysical models of massive and stellar-mass black hole populations to compute, for the first time, the stochastic gravitational-wave background from these sources. In optimistic scenarios the background is observable by Advanced LIGO and LISA for field masses $m_s$ in the range $\sim [2\times 10^{-13}, 10^{-12}]\,{\rm eV}$ and $\sim 5\times[ 10^{-19}, 10^{-16}]\,{\rm eV}$, respectively, and it can affect the detectability of resolvable sources. Our estimates suggest that an analysis of the stochastic background limits from LIGO O1 might already be used to marginally exclude axions with mass $\sim 10^{-12.5}{\rm eV}$. Semicoherent searches with Advanced LIGO (LISA) should detect $\sim 15~(5)$ to $200~(40)$ resolvable sources for scalar field masses $3\times 10^{-13}$ ($10^{-17}$) eV. LISA measurements of massive BH spins could either rule out bosons in the range $\sim [10^{-18}, 2\times 10^{-13}]$ eV, or measure $m_s$ with ten percent accuracy in the range $\sim[10^{-17}, 10^{-13}]$ eV.

Equations

$${\tau }_{\mathrm{inst}}$$

$${\tau }_{\mathrm{GW}}$$

$$h=\sqrt{\frac{2}{5\pi }}\frac{GM}{c^{2}r}\left(\frac{M_S}{M}\right)A(\chi ,f_{s}M),$$

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

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

$$\begin{gathered} {\Omega }_{\mathrm{GW}}(f)=\frac{f}{{\rho }_c}{\int }_{\rho <8}dz\frac{dt}{dz} \\ \dot{n}(M,\chi ,z)\frac{d{E}_s}{d{f}_s}, \end{gathered}$$

$${\rho }_{\mathrm{stoch}}=\sqrt{T_{\mathrm{obs}}{\int }_{f_{\min }}^{f_{\max }}df\frac{{\Omega }_{\mathrm{GW}}^2}{{\Omega }_{\mathrm{sens}}^2}},$$

$$N={\int }_{\rho >8}\frac{d^2\dot{n}}{dMd\chi }\left(\frac{T_{\mathrm{obs}}}{1+z}+\Delta t\right)\frac{d{V}_c}{dz}dzdMd\chi ,$$

Full text

Introduction. The historical LIGO gravitational wave (GW) detections Abbott et al. (2016a, b, 2017a) provide the strongest evidence to date that astrophysical black holes (BHs) exist and merge Berti et al. (2006); Cardoso et al. (2016); Berti et al. (2016). Besides probing the nature of compact objects and testing general relativity Gair et al. (2013); Yunes and Siemens (2013); Berti et al. (2015); Abbott et al. (2016c), LIGO Aasi et al. (2015) and the space-based detector LISA Audley et al. (2017) may revolutionize our understanding of particle physics and dark matter. Ultralight bosons, which could be a significant component of dark matter Arvanitaki et al. (2010); Essig et al. (2013); Marsh (2016); Hui et al. (2016), interact very weakly (if at all) with baryonic matter, but the equivalence principle implies that their gravitational interaction should be universal. Low-energy bosons near spinning BHs can trigger a superradiant instability whenever the boson frequency $\omega_{R}$ satisfies the superradiant condition $0<\omega_{R}<m\Omega_{\rm H}$, where $\Omega_{\rm H}$ is the horizon angular velocity and $m$ is an azimuthal quantum number, with possible astrophysical implications Press and Teukolsky (1972); Detweiler (1980); Cardoso et al. (2004); Brito et al. (2015a).

Despite extensive work on massive spin-0 Detweiler (1980); Zouros and Eardley (1979); Cardoso and Yoshida (2005); Dolan (2007), spin-1 Pani et al. (2012a, b); Witek et al. (2013); Endlich and Penco (2016); East and Pretorius (2017); Baryakhtar et al. (2017) and spin-2 fields Brito et al. (2013), the evolution and the end-state of the instability are not fully understood Okawa et al. (2014); Zilhão et al. (2015); Sanchis-Gual et al. (2016); Bosch et al. (2016). Recent numerical simulations East and Pretorius (2017) support the conclusions of perturbative studies Arvanitaki and Dubovsky (2011); Brito et al. (2015b); Arvanitaki et al. (2015); Yoshino and Kodama (2014a); Brito et al. (2015a); Arvanitaki et al. (2017); Baryakhtar et al. (2017): the BH spins down, transferring energy and angular momentum to a mostly dipolar boson condensate until $\omega_{R}\sim m\Omega_{H}$. The energy scale is set by the boson mass $m_{s}\equiv\mu\hbar$, which implies that $\omega_{R}\sim\mu$ and that the instability saturates at $\mu\sim m\Omega_{\rm H}$ (in units $G=c=1$). The condensate is then dissipated through the emission of mostly quadrupolar GWs, with frequency set by $\mu$. The mechanism is most effective when the boson’s Compton wavelength is comparable to the BH’s gravitational radius: detailed calculations show that the maximum instability rate for scalar fields corresponds to $M\mu\simeq 0.42$ Dolan (2007). Therefore, the instability window corresponds to masses $m_{s}\sim 10^{-14}$–$10^{-10}$ eV and $m_{s}\sim 10^{-19}$–$10^{-15}$ eV for LIGO and LISA BH-boson condensate sources, respectively Brito et al. (2015a). In this work and in a companion paper Brito et al. (2017) we argue that GW detectors can discover new particles beyond the Standard Model or impose constraints on their masses.

GWs from scalar condensates around BHs. The instability occurs in two stages Brito et al. (2015b). In the first (linear) phase the condensate grows on a timescale $\tau_{\rm inst}\sim M^{-8}\mu^{-9}$ until the superradiant condition is nearly saturated. In the second (nonlinear) phase GW emission governs the evolution of the condensate, which is dissipated over a timescale $\tau_{\rm GW}$ that depends on its mass $M_{S}$ and on the GW emission rate. These two timescales can be computed analytically when $M\mu\ll 1$ Brito et al. (2017). For small dimensionless BH spins $\chi\equiv J/M^{2}\ll 1$, they read

[TABLE]

These relations (valid for any BH mass) are a good approximation even when $M\mu$ and $\chi$ are $\sim 1$ Brito et al. (2017). Since $\tau_{\rm GW}\gg\tau_{\rm inst}\gg M$, the condensate has enough time to grow, and the evolution of the system can be studied in a quasi-adiabatic approximation Brito et al. (2015b) using Teukolsky’s formalism Teukolsky (1973); Yoshino and Kodama (2014b). The field’s stress-energy tensor is typically small, thus its backreaction is negligible Brito et al. (2015b); East and Pretorius (2017).

Over the emission timescale (which in most cases is much longer that the observation time $T_{\rm obs}$), the GWs are nearly monochromatic, with frequency $f_{s}=\omega_{R}/\pi\sim\mu/\pi$. As such, BH-boson condensates are continuous sources, like pulsars for LIGO or verification binaries for LISA. We conservatively assume that GWs are produced after saturation of the instability, which leads the BH from an initial state $(M_{i},\,J_{i})$ to a final state $(M,\,J)$, and we thus compute the root-mean-square strain amplitude $h$ using the final BH parameters. By averaging over source and detector orientations we get

[TABLE]

where $r$ is the (comoving) distance to the source, the masses are in the source frame, and the dimensionless function $A(\chi,f_{s}M)$ is computed from BH perturbation theory Yoshino and Kodama (2014b); Brito et al. (2017). Our results are more accurate than the analytic approximations of Arvanitaki and Dubovsky (2011); Brito et al. (2015b). It can be shown that $M_{S}$ scales linearly with $J_{i}$ Brito et al. (2017), so $h$ also grows with $J_{i}$. For LISA, we also take into account correction factors due to the detector geometry Berti et al. (2005). In the detector frame, Eq. (3) still holds if the masses $M$ and $M_{S}$ are multiplied by $(1+z)$, $r$ is replaced by the luminosity distance, and the frequency is replaced by the detector-frame frequency $f=f_{s}/(1+z)$. Nevertheless, one needs to use detector-frame frequencies when comparing to the detector sensitivity.

In semicoherent searches of monochromatic sources, the signal is divided in ${\cal N}$ coherent segments of time length $T_{\rm coh}$, and we have $h_{\rm thr}\simeq 25{\cal N}^{-1/4}\sqrt{S_{h}(f)/T_{\rm coh}}$, where $h_{\rm thr}$ is the minimum root-mean-square strain amplitude detectable over the observation time ${\cal N}\times T_{\rm coh}$ Palomba (2012), and $S_{h}(f)$ is the noise power spectral density (PSD) at $f$ Ruiter et al. (2010).

In Fig. 1 we compare the GW strain of Eq. (3) with the PSDs of LISA and Advanced LIGO at design sensitivity. The GW strain increases almost vertically as a function of $\omega_{R}\simeq\mu$ in the superradiant range $(0,\Omega_{H})$. Thin solid curves correspond to the stochastic background from the whole BH population, for a boson mass $m_{s}$. This background produces itself a “confusion noise” when $m_{s}\approx[10^{-18},10^{-16}]$ eV, complicating the detection of individual sources. Figure 1 suggests that bosons with masses $10^{-19}\,{\rm eV}\lesssim m_{s}\lesssim 10^{-11}\,{\rm eV}$ (with a small gap around $m_{s}\sim 10^{-14}$ eV, which might be filled by DECIGO Kawamura et al. (2006)) could be detectable by LIGO and LISA. Below we quantify this expectation.

BH population models. Assessing the detectability of these signals requires astrophysical models for BH populations. For LISA sources, the main uncertainties concern the mass and spin distribution of isolated BHs, the model for their high-redshift seeds, and their accretion and merger history. We adopt the same populations of Klein et al. (2016); Babak et al. (2017), which were based on the semianalytic galaxy formation calculations of Barausse (2012) (see also Sesana et al. (2014); Antonini et al. (2015a, b)). In our optimistic model, we use these calculations to infer the redshift-dependent BH number density ${d^{2}n}/(d\log_{10}Md\chi)$. The spin distribution is skewed toward $\chi_{i}\sim 1$, at least at low masses Sesana et al. (2014). We also adopt less optimistic and pessimistic models with mass function given by Eqs. (5) and (6) of Babak et al. (2017) for $z<3$ and $10^{4}M_{\odot}<M<10^{7}M_{\odot}$, whereas for $M>10^{7}M_{\odot}$ we use a mass distribution with normalization $10$ and $100$ times lower than the optimistic one. In both the less optimistic and pessimistic models we assume a uniform spin distribution in the range $\chi_{i}\in[0,1]$.

The LIGO stochastic GW background comes mostly from extra-galactic stellar-mass BHs, which were ignored in previous work Arvanitaki et al. (2015). Here we model these sources using the semianalytic galaxy evolution model of Dvorkin et al. (2016). The BH formation rate as a function of mass and redshift reads

[TABLE]

where $\tau({\cal M_{\star}})$ is the lifetime of a star of mass ${\cal M_{\star}}$, $\phi({\cal M_{\star}})$ is the stellar initial mass function, $\psi(t)$ is the cosmic star formation rate (SFR) density and $\delta$ is the Dirac delta. We fit the cosmic SFR as described in Vangioni et al. (2015) and calibrate it to observations of luminous galaxies Bouwens et al. (2011); Behroozi et al. (2013a). We assume a Salpeter initial mass function $\phi({\cal M_{\star}})\propto{\cal M_{\star}}^{-2.35}$ Salpeter (1955) in the range ${\cal M_{\star}}\in[0.1-100]\,M_{\odot}$, and take stellar lifetimes from Schaerer (2002). We also follow the production of metals by stars Woosley and Weaver (1995) and the resulting enrichment of the interstellar medium, which affects the metallicity of subsequent stellar generations. The function $g({\cal M_{\star}})$ relates the initial stellar mass $\cal M_{\star}$ and the BH mass $M$, and encodes the BH formation process. In general, the mass of the BH formed from a star with initial mass ${\cal M_{\star}}$ depends on the stellar metallicity (Vink, 2008) and rotational velocity (de Mink et al., 2009), as well as interactions with its companion if the star belongs to a binary system. We assume that all stellar-mass BHs are produced from isolated massive stars after core collapse, and calculate the BH mass for a given ${\cal M_{\star}}$ and metallicity using the analytic fits for the “delayed” model of Fryer et al. (2012). Through the metallicity, the function $M=g({\cal M_{\star}})$ is implicitly a function of redshift. Since this model does not predict the initial BH spins, we assume a uniform distribution and explore different ranges: $\chi_{i}\in[0.8,1]$, $[0.5,1]$, $[0,1]$ and $[0,0.5]$.

The dominant contribution to LIGO resolvable signals comes from Galactic stellar-mass BHs Arvanitaki et al. (2015). We estimate their present-day mass function as

[TABLE]

where the integration is over all cosmic times prior to the present epoch; $N_{\rm MW}$ denotes the number of BHs in the Galaxy; ${\rm SFR}(z)$ is the SFR of Milky-Way type galaxies as a function of redshift Behroozi et al. (2013a, b); ${dp}/{d{\cal M_{\star}}}$ is the probability of forming a star with mass between ${\cal M_{\star}}$ and ${\cal M_{\star}}+d{\cal M_{\star}}$ (obtained from the Salpeter initial mass function); and $dM/{d{\cal M_{\star}}}$ is given by the “delayed” model of Fryer et al. (2012). This latter quantity is a function of redshift through the metallicity, whose redshift evolution we model following Ma et al. (2016). To obtain a (differential) BH number density ${dn_{\rm MW}}/{dM}$, we “spread” this mass function over the Galaxy, proportionally to the (present) stellar density. For the latter we assume a simple bulge+disk model, where the bulge is modeled via a Hernquist profile Hernquist (1990) with mass $\sim 2\times 10^{10}M_{\odot}$ and scale radius $\sim 1\,{\rm kpc}$ Shen et al. (2003), and the disk has an exponential profile with mass $\sim 6\times 10^{10}M_{\odot}$ and scale radius $\sim 2$ kpc Porcel et al. (1998).

Stochastic background. The stochastic background produced by BH-boson condensates is given by an integral over unresolved sources – those with signal-to-noise ratio (SNR) $\rho<8$ – of the formation rate density per comoving volume $\dot{n}$ Phinney (2001):

[TABLE]

where $\rho_{c}=3H_{0}^{2}/(8\pi)\approx 1.3\times 10^{11}M_{\odot}/{\rm Mpc}^{3}$ is the critical density of the Universe, ${{\rm}dt}/{{\rm}dz}$ is the derivative of the lookback time $t(z)$ with respect to $z$, ${{\rm}dE_{s}}/{{d}f_{s}}$ is the energy spectrum in the source frame, and $f$ is the detector-frame frequency. For LIGO we compute $\dot{n}$ by integrating Eqs. (4) and (5). For LISA we integrate ${d^{2}n}/(d\log_{10}Md\chi)$ – as given by the aforementioned “optimistic”, “less optimistic” and “pessimistic” models – with respect to mass and spin, and we assume that $\dot{n}=n/t_{0}$, where $t_{0}\approx 13.8\,{\rm Gyr}$ is the age of the Universe (i.e., each BH undergoes boson annihilation only once in its cosmic history). This assumption does not significantly affect our results, because subsequent annihilation signals (if they occur at all) are much weaker Brito et al. (2017).

For the spectrum of the GW signal we assume ${{\rm}dE_{s}}/{{d}f_{s}}\approx E_{\rm GW}\delta(f(1+z)-f_{s})$, where $E_{\rm GW}$ is the total energy radiated in GWs over the signal duration $\Delta t$, and the Dirac delta is “spread out” over a frequency window of width $\sim\max[1/(\Delta t(1+z)),1/T_{\rm obs}]$ to account for the finite signal duration and observation time. For LIGO we can safely neglect the effect of mergers Gerosa and Berti (2017); Fishbach et al. (2017) and accretion King and Kolb (1999). For LISA, we conservatively assume that mergers and accretion cut the signal short, and thus define the signal duration as $\Delta t=\left\langle\min\left(\tau_{\rm GW}/(N_{m}+1),\,t_{S},\,t_{0}\right)\right\rangle$, where $\tau_{\rm GW}$ is given by Eq. (2); $t_{S}=4.5\times 10^{8}{\rm\,yr\,}{\eta}/[{f_{\rm Edd}(1-\eta)}]$ is the typical accretion “Salpeter” timescale, which depends on the Eddington ratio $f_{\rm Edd}$ and on the spin-dependent radiative efficiency $\eta$; $\langle...\rangle$ denotes an average weighted by the Eddington-ratio probability distribution; and $N_{m}$ is the average number of mergers in the interval $[t(z)-\frac{1}{2}\tau_{\rm GW},t(z)+\frac{1}{2}\tau_{\rm GW}]$ Brito et al. (2017). Moreover, since our calculation assumes that the instability saturates before GW emission takes place, our stochastic background calculation only includes BHs for which the expected number of mergers during the instability timescale is $N_{m}<1$, and for which $\tau_{\rm inst}<\Delta t$ (thus ensuring that the instability timescale is shorter than the typical accretion and merger timescales).

The ${\rm SNR}$ for the stochastic background is Thrane and Romano (2013)

[TABLE]

where $\Omega^{\rm LIGO}_{\rm sens}=\frac{S_{h}(f)}{\sqrt{2}\Gamma_{IJ}(f)}\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}$ and $\Omega^{\rm LISA}_{\rm sens}=S_{h}(f)\frac{2\pi^{2}}{3H_{0}^{2}}f^{3}$ for LIGO Allen and Romano (1999) and LISA Cornish and Larson (2001), respectively. In the LIGO case we assume the same $S_{h}$ for the Livingston and Hanford detectors, and $\Gamma_{IJ}$ denotes their overlap reduction function Thrane and Romano (2013).

The order of magnitude of the stochastic background shown in Fig. 2 (left panel) can be estimated by a simple back-of-the-envelope calculation. The average mass fraction of an isolated BH emitted by the boson cloud is $f_{\rm ax}\sim{\cal O}(1\%)$ Brito et al. (2017). Because the signal is almost monochromatic, the emitted GWs in the detector frame span about a decade in frequency, i.e. $\Delta\ln f\sim 1$ for both LISA and LIGO (cf. Fig. 2). Thus, $\Omega_{\rm GW,\,ax}=({1}/{\rho_{\rm c}})({d\rho_{\rm GW}}/{d\ln f})\sim f_{\rm ax}{\rho_{\rm BH}}/{\rho_{\rm c}}$, where $\rho_{\rm GW}$ and $\rho_{\rm BH}$ are the GW and BH energy density, respectively. Since the BH mass density is $\rho_{\rm BH}\sim{\cal O}(10^{4})M_{\odot}/{\rm Mpc}^{3}$ in the mass range $10^{4}-10^{7}M_{\odot}$ relevant for LISA, this yields $\Omega^{\rm LISA}_{\rm GW,\,ax}\sim 10^{-9}$. For LIGO, the background of GWs from BH binaries can be approximated as $\Omega_{\rm GW,\,bin}\sim f_{\rm GW}f_{\rm m}\rho_{\rm BH}/\rho_{c}$, where $f_{\rm GW}\sim{\cal O}(1\%)$ is the binary’s mass fraction emitted in GWs Barausse et al. (2012), and $f_{\rm m}\sim{\cal O}(1\%)$ Dvorkin et al. (2016) is the fraction of stellar-mass BHs in binaries that merge in less than $t_{0}$. Therefore $\Omega_{\rm GW,\,ax}/\Omega_{\rm GW,\,bin}\sim f_{\rm ax}/(f_{\rm GW}f_{\rm m})\sim 10^{2}$. Since the O1 results imply peak background values $\Omega_{\rm GW,\,bin}\sim 10^{-9}-10^{-8}$ Abbott et al. (2016d, 2017b) (or larger if spins are included), we obtain $\Omega^{\rm LIGO}_{\rm GW,\,ax}\sim 10^{-7}-10^{-6}$. These estimates are in qualitative agreement with the left panel of Fig. 2.

Remarkably, $\rho_{\rm stoch}$ (right panel of Fig. 2) can be very high. For optimistic astrophysical models, boson masses in the range $2\times 10^{-13}\,{\rm eV}\lesssim m_{s}\lesssim 10^{-12}\,{\rm eV}$ ($5\times 10^{-19}\,{\rm eV}\lesssim m_{s}\lesssim 5\times 10^{-16}\,{\rm eV}$) yield $\rho_{\rm stoch}>8$ with LIGO (LISA). Our estimates suggest that, for the most pessimistic model and masses around $m_{s}\approx 3\times 10^{-12}\,{\rm eV}$, the background would have ${\rm SNR}\approx 1.2$ using our simple analytic estimate of the LIGO O1 sensitivity, thus being only marginally allowed by current LIGO O1 upper limits Abbott et al. (2017b). Our conclusions should be validated by a careful data analysis of the stochastic background in LIGO O1 and O2. In particular, current upper limits on the stochastic background assume that the spectrum can be described by a power law in the LIGO range Abbott et al. (2017b), which is not the case for the backgrounds computed here.

Resolvable sources. We estimate the number of resolvable events as Brito et al. (2017)

[TABLE]

where $dV_{c}=4\pi D_{c}^{2}{dD_{c}}$, $D_{c}$ is the comoving distance, and $\dot{n}=n/t_{0}$ for LISA. The dependence on $T_{\rm obs}/(1+z)+\Delta t$ comes about because the probability that an observation of duration $T_{\rm obs}$ and a signal of duration $\Delta t(1+z)$ (in the detector frame) overlap is proportional to the sum of the two durations. In the limit $\Delta t(1+z)\ll T_{\rm obs}$ we have $N\propto T_{\rm obs}$, as usual for short-lived sources Hartwig et al. (2016). For $\Delta t(1+z)\gg T_{\rm obs}$, $N$ becomes proportional to the duty cycle $\Delta t/t_{f}$, $t_{f}\equiv n/\dot{n}$ being the formation timescale of the boson condensates. This duty cycle, akin e.g. to the duty cycle of active galactic nuclei, accounts for the fact that only a fraction of the sources are radiating during the observation time.

Figure 3 shows resolvable event rates assuming (conservatively) semicoherent searches for different astrophysical models. The SNR was computed by including the confusion noise from the stochastic background of unresolvable boson-condensate sources: neglecting this contribution would overestimate LISA rates by more than one order of magnitude. Our models typically predict $\sim$ 40 (200) events in $121\times 250$ hours of total observation time for the optimistic models and boson masses in the optimal range around $m_{s}\sim 10^{-17}\,{\rm eV}$ ($3\times 10^{-13}\,{\rm eV}$) for LISA (LIGO). Rates in the less optimistic and pessimistic models decrease by factors of order unity. However, it is remarkable that a boson with $m_{s}\sim 10^{-17}\,{\rm eV}$ ($3\times 10^{-13}\,{\rm eV}$) would produce around $5$ ($15$) direct LISA (LIGO) detections even for pessimistic astrophysical models.

So far we focused on the direct detection of GWs from bosonic condensates. In Brito et al. (2017) we use Bayesian model selection to show that LISA could infer the existence of light bosons indirectly: LISA measurements of massive BH spins could provide evidence for holes in the BH mass-spin “Regge plane” (i.e., for the absence of BHs spinning above the superradiant instability window) Arvanitaki and Dubovsky (2011). As indicated by the shaded areas in Fig. 3, a 4-year LISA mission could rule out boson masses in a range that depends on the assumed BH model ($[4.5\times 10^{-18},1.6\times 10^{-13}]$ for the “light-seed” popIII model, $[10^{-18},2.3\times 10^{-14}]$ for the “heavy-seed, no-delay” Q3-nod model of Klein et al. (2016)). If fields with $m_{s}\in[10^{-17},10^{-13}]$ eV exist in nature, LISA observations of BH mergers can measure $m_{s}$ with ten percent accuracy Brito et al. (2017).

Conclusions. Together, Earth- and space-based detectors will allow for multiband GW searches of ultralight bosons in the range $[10^{-19}-10^{-10}]\,{\rm eV}$. We plan to improve estimates of the stochastic background for LIGO by using population synthesis models Berti et al. (2016); Belczynski et al. (2016a, b). The potential of detectors like DECIGO or the Einstein Telescope to detect or rule out bosons of mass $m_{s}\sim 10^{-14}\,{\rm eV}$ should also be investigated by using intermediate-mass BH formation models Kawamura et al. (2006); Gair et al. (2011). Our analysis must be extended to spin-1 Pani et al. (2012a, b) and spin-2 Brito et al. (2013) fields, for which the instability time scales are shorter and GW amplitudes are larger. Our results also suggest that recent estimates of resolvable GWs from spin-1 instabilities Baryakhtar et al. (2017) should be revised taking into account the stochastic background.

Acknowledgments. We would like to thank Lijing Shao for useful comments. S. Ghosh and E. Berti are supported by NSF Grants No. PHY-1607130 and AST-1716715. E. Berti is supported 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. 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. 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. 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).