The stochastic gravitational-wave background in the absence of horizons

According to general relativity and the standard model of particle physics, dark and compact objects more massive than  ≈$3M_{\odot}$ must be black holes (BHs) [1]. These are characterized by an event horizon, causally disconnecting the BH interior from its exterior, where observations take place. In classical gravity, BHs can form from gravitational collapse [2], providing a sound and compelling theoretical support for their existence. When quantum effects are included at semiclassical level, however, BHs are not completely 'dark', but evaporate by emitting thermal black-body radiation [3]. For astrophysical BHs, evaporation is negligible. Therefore, BHs are commonly accepted to exist—with masses in the ranges from  ∼$10 M_\odot$ to  ∼$60 M_\odot$ (and perhaps larger mass), and from  ∼$10^{5} M_\odot$ to  ∼$10^{10} M_\odot$—and play a fundamental role in astronomy and astrophysics [4].

It is sometimes not fully appreciated that BHs are truly 'holes' in spacetime, where time 'ends' and inside which the known laws of classical physics break down [5]. Furthermore, the classical concept of an event horizon seems at clash with quantum mechanics, and the very existence of BHs leads to unsolved conundra such as information loss [5]. Thus, in reality, the existence of BHs is an outstanding event, for which one should provide equally impressive evidence [6–8]. In previous decades several alternatives and arguments have been put forward according to which—in a quantum theory—BHs would either not form at all, or would just be an ensemble of horizonless quantum states [9–15].

From a theoretical standpoint, BHs therefore lay at the interface between classical gravity, quantum theory and thermodynamics, and understanding their nature may provide a portal to quantum gravity or other surprises. In the formal mathematical sense, it is impossible to ever show that BHs exist, since in general relativity their definition requires knowledge of the whole spacetime, including the future [2]. However, the newborn gravitational-wave (GW) astronomy allows us to constrain alternatives to BHs to unprecedented level. GW detectors can rule out a wide range of models, through observations of inspiralling binaries or the relaxation of the final object forming from a merger [6–8, 16–18]. Here, we explore one significant effect that follows from the absence of the most salient feature of a BH, the event horizon. We will show that compact, horizonless spinning geometries would fill the universe with a background of GWs detectable by current and future instruments, through a classical process known as the 'ergoregion instability' [19, 20].

In Einstein's theory, the unique globally vacuum astrophysical solution for a spinning object is the Kerr geometry. It depends on two parameters only: its mass M and angular momentum $J=GM^2\chi/c$, with G Newton's constant, c the speed of light and $|\chi|\leqslant1$ a pure number. The compact, dark objects in our universe could depart from the Kerr geometry in two distinct ways. The near-horizon structure might change significantly, while retaining the horizon [21–23]. In coalescing binaries, such effects can be probed by GW measurements of the quadrupole moment, the tidal absorption and deformability [7, 8, 24, 25], and especially the quasinormal oscillation modes of the remnant object [26, 27]. Here we explore a second (and more subtle) scenario, where the geometry is nearly everywhere the same as that of a BH, but the horizon is absent. Two smoking-gun effects arise in this scenario. First, the late-time ringdown consists generically of a series of slowly damped 'echoes' [6, 8, 25, 28]. Furthermore, by working as a one-way membrane, horizons act as a sink for external fluctuations, including those inside the ergoregion, where negative-energy states are possible [20, 29]. Such states are typically associated with instabilities: their existence allows scattering waves to be amplified, i.e. positive-energy perturbations can be produced, which can travel out of the ergoregion. Energy conservation then requires the negative-energy states inside the ergoregion to grow. In a BH, this piling up can be avoided by dumping the negative energy into the horizon, thus stabilizing the object. In the absence of a horizon, instead, this process leads to an exponential cascade. As a consequence, spinning BHs are linearly stable, but any horizonless object sufficiently similar to a rotating BH is unstable [19, 30–32].

2.1. Canonical model: perfectly reflecting surface

We start by the simplest model of horizonless geometries: a compact body whose exterior is described by the Kerr metric, and with a perfectly reflective surface. This spacetime defines a natural cavity, i.e. the region between the object's surface and the potential barrier for massless particles (the 'photon sphere') (see figure 1). In this cavity, negative-energy modes, and thus instabilities, can be excited [15, 33]. The dynamics is controlled by two parameters: the size of the cavity and the object's angular velocity, which determine how fast the instability grows. The cavity size can be parametrized by the light travel time t₀ (as observed at infinity) between the photon sphere and the object's surface [6–8, 20]. The timescale t₀ also defines a set of possible modes, with fundamental frequencies $\omega=\omega_R+ {\rm i}\omega_I$ and $\omega_R\sim \pi/t_0$. The instability is controlled by the amplification factor $|{{\mathcal A}}|^2$ of the ergoregion at this frequency [20], i.e. $\omega_I\sim |{{\mathcal A}}|^2/t_0$. This follows from a very generic 'bounce-and-amplify' argument, which was shown to accurately describe specific models [15, 20].

Zoom In Zoom Out Reset image size

A scalar, electromagnetic, or gravitational perturbation in such a geometry grows exponentially on a timescale $\newcommand{\e}{{\rm e}} \tau\equiv 1/\omega_I$. The characteristic unstable modes can be computed numerically and agree well with bounce-and-amplify estimates [15, 34]. When $\newcommand{\e}{{\rm e}} t_0\gg t_H\equiv {GM}/{c^3}$, these modes are well described by [15, 32–35]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_R \simeq -\frac{q\pi}{2t_0}+m\Omega , \label{omegaR} \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_I \simeq -\beta_{ls}\frac{t_H}{t_0}(1+\Delta){\Delta^{2l}} \left(2\omega_R t_H\right)^{2l+1} (\omega_R-m\Omega), \label{omegaI} \nonumber \end{align} \tag{ 2 }$$

where $\Delta=\sqrt{1-\chi^2}$, $\Omega(\chi)$ is the object's angular velocity, ${\beta_{ls}^{1/2}}={(l-s)!(l+s)!}/$$[{(2l)!(2l+1)!!}]$ for a spin-s perturbation with angular number l [20, 36], and q is a positive odd (even) integer for Dirichlet (Neumann) boundary conditions at the surface [34]. Thus, the spacetime is unstable for $\omega_R(\omega_R-m\Omega)<0$ (i.e. in the superradiant regime [20]), on a timescale $\newcommand{\e}{{\rm e}} \tau\equiv 1/\omega_I$. Here we consider the dominant gravitational mode (l  =  m  =  2) to estimate the background of GWs.

Equations (1) and (2) are valid for any object able to completely reflect the incoming radiation. In particular, if the object's surface sits at a constant (Boyer–Lindquist) radius

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_0=r_+(1+\epsilon), \nonumber \end{align} \tag{ 3 }$$

where $r_+=GM(1+\Delta)/c^2$ is the location of the (would-be) event horizon in these coordinates, then the travel time reads [6, 8, 25]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle t_0^{\rm canonical} \sim t_H |\log\epsilon| .\label{t0_canonical} \nonumber \end{align} \tag{ 4 }$$

Several different arguments about the magnitude of can be made. If quantum-gravity effects become important at Planck timescales $t_{\rm P}=\sqrt{\hbar G/c^5}$, it is natural to set $\newcommand{\e}{{\rm e}} \epsilon=t_{\rm P}/t_H\sim 10^{-39}{{\rm \mbox{--}}}10^{-46}$ for stellar-mass to supermassive dark objects. These objects were dubbed ClePhOs in [15], and are impossible to rule out in practice via electromagnetic observations [8, 15].

2.2. Modelling the interior

2.3. Evolution of the instability

The evolution of the object's mass and angular momentum, under energy and angular momentum losses, can be computed within the adiabatic approximation (because $\tau\gg t_H$) [39]. The unstable mode is simply draining energy and angular momentum from the object, which we assume to have an equation of state such that $\newcommand{\e}{{\rm e}} \epsilon={\rm const}$ during the evolution. From energy and angular momentum conservation, the evolution equations for each mode read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{evol_equations} \dot E = \dot{E}_0\,e^{2\omega_I t}, \quad \dot{M}=-\dot E/c^2,\quad \dot{J}=-m{\dot E}/{\omega_R}, \nonumber \end{align} \tag{ 5 }$$

where $\dot{E}_0$ encodes the initial preturbation of the (unstable) system. Since the instability is exponential, the overall evolution is insensitive to the precise value of these initial conditions. The equations above are valid for a monochromatic mode in a generic stationary and axisymmetric background. The energy flux can be written as ${\rm d}E/{\rm d}f = \dot E /\dot f$, where $f=\omega_R/(2\pi)$ is the frequency associated with the mode. From the evolution equations, we can evaluate $\dot M$ and $\dot J$ and, in turn, $\dot \omega_R$. To leading order in the $\newcommand{\e}{{\rm e}} \epsilon\to0$ limit, we obtain $\omega_R\sim m\Omega$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{dEdf} \frac{{\rm d}E}{{\rm d}f} \sim \frac{16G^2\pi^2}{m^2c^4} f M^3, \nonumber \end{align} \tag{ 6 }$$

valid for any angular numbers l and m. In the same limit, the critical value of spin above which the ergoregion instability occurs reads [34, 35]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_{\rm crit} \approx\frac{\pi q}{m |\log(\epsilon)|} \approx \frac{0.035\,q}{m} \log^{-1}\left(\frac{\epsilon}{10^{-40}}\right). \nonumber \end{align} \tag{ 7 }$$

Thus, if $\chi(t=0)>\chi_{\rm crit}$, the instability removes energy and angular momentum until superradiance is saturated, i.e. $\chi(t\gg\tau)\to\chi_{\rm crit}$. Note that the small spin value $\chi_{\rm crit}$ is compatible with the low measured spins of inspiralling compact objects detected via GWs so far [40]. Since we are interested in gravitational perturbations, when solving the evolution equations (5) and computing the energy flux (6) we only consider the dominant l  =  m  =  2 mode and neglect higher-modes. For each initial spin χ and travel time t₀ we only consider the dominant q-mode associated with the shortest instability timescale. Note also that the above analysis assumes that the backreaction of these fields on the geometry is negligible. Our results indicate that those are always reliable approximations.

2.4. GW stochastic background

A population of GW sources too far and/or weak to be detected individually may still give rise to a 'stochastic' background detectable by a network of interferometers, e.g. the LIGO/Virgo network, sensitive to frequency ranging from  ∼10 to  ∼100 Hz [41, 42]; the pulsar-timing-array experiments [43–47], which are already constraining backgrounds at frequencies  ∼10⁻⁹–10⁻⁶ Hz; and the future LISA constellation [48], which will be sensitive to frequencies between 10⁻⁶ Hz and  ∼1 Hz. The background is produced by the incoherent superposition, at the detector, of the GW signals from all the unresolved sources in the population. The background can be characterized either by (i) its (dimensionless) energy spectrum

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{\rm gw}(\,f_o)=\frac{1}{\rho_{\rm c}}\frac{{\rm d}\rho_{\rm gw}}{{\rm d}\ln f_o}, \label{OmegaGW} \nonumber \end{align} \tag{ 8 }$$

($\rho_{\rm gw}$ being the background's energy density, f_o the frequency measured at the detector and $\rho_{\rm c}$ the critical density of the Universe at the present time), obtained by summing the energies emitted by all the unresolved sources in a given frequency bin [49]; or (ii) directly by the characteristic strain $h_{\rm c}(\,f_o)$ observed in the detector, which can be obtained summing in quadrature (and binning in frequency) the strain amplitudes of all the unresolved sources [50]. The two quantities are related by $\Omega_{\rm gw}(\,f_o)={2\pi^2} \left[\,f_o h_{\rm c}(\,f_o)\right]^2/{(3H_0^2)}$, where $H_0\approx 68~{\rm km~(s~Mpc){}^{-1}}$ is the Hubble rate. While these two ways of computing the background signal are equivalent [50], we have implemented both as a consistency check of our results. This also allows us to check that the number of sources contributing in each frequency bin is typically large as long as the bin size is $\gtrsim$0.01 dex in the LISA band (which ensures that the number of sources contributing 99% of the signal in each bin ranges from thousands to millions). In the LIGO band, sources are even more numerous: frequency bins $\gtrsim$0.01 dex yield 10⁹–10¹⁴ sources contributing 99% of the signal in each bin. (These are mostly extragalactic sources, as Galactic ones give a negligible contribution to the background.) This in turn implies, in particular, that the background is expected to be smooth with that frequency binning [50]. When computing the stochastic background of unstable exotic compact objects, the energy flux given by equation (6) is defined in the frequency range $f\in [\,f_{\rm min}, f_{\rm max}]$; f_max can be computed using equation (1) for a given initial mass and spin of the compact object, and f_min is computed by solving the evolution equation (5) from the formation redshift of the compact object to the present time.

For the astrophysical populations of isolated BHs, we adopt the same models as in [49, 54]. We assume that all such astrophysical BH candidates are actually BH mimickers. For stellar-origin BHs we account for both Galactic and extragalactic BHs that form from the core collapse of massive ($\gtrsim$$20 M_\odot$) stars, by tracking the cosmic star formation history and the metallicity evolution of the Universe [55]. We assume a uniform distribution for the initial spins with $\chi\in[0, 1]$ as the most optimistic and $\chi\in[0, 0.5]$ as the most pessimistic scenario. For the massive (∼$10^4{{\rm \mbox{--}}}10^7 M_\odot$) and supermassive (∼$10^8{{\rm \mbox{--}}}10^{10} M_\odot$) BHs that emit respectively in the LISA and PTA bands, we adopt the semi-analytic galaxy-formation model of [56] (with later incremental improvements described in [57–59]), which follows the formation of these objects from their high-redshift seeds and their growth by accretion and mergers. This growth is triggered in turn by the synergic co-evolution of the BHs with their host galaxies, of which we evolve both the various baryonic components and the dark-matter halos. This model is optimistic since it predicts a spin distribution skewed towards large spins, at least at low masses. To include astrophysical uncertainties in our computation, we also consider models in between our most optimistic and most pessimistic assumptions as described in [49] (see section 3 therein).

Our main results for the GW stochastic background from exotic compact objects are shown in figure 2 in the frequency bands relevant for LIGO/Virgo (top panel) and for LISA/an SKA-based pulsar timing arrays (bottom panel). The top panel suggests that the absence of a stochastic background in LIGO O1 already rules out our canonical model even for conservative spin distributions, while LIGO at design sensitivity will be able to rule out our canonical model even in more pessimistic scenarios than those assumed here, e.g. even if all BH-like objects had initial spin $\chi<0.2$. Similar results apply in the LISA band, whereas the stochastic signal is too small to be detectable by pulsar timing arrays, even in the SKA era⁶.

Zoom In Zoom Out Reset image size

The level of the stochastic background shown in figure 2 can also be understood with an approximate analytic calculation [54]. The BH-mimicker mass fraction lost to GWs due to superradiance is $F_{\rm sr}\sim {{\mathcal O}}(1\%)$ [49, 61]. Because the signal spans about a decade in frequency (see equation (1)), $\Delta \ln f \sim 1$, and $\Omega_{\rm GW, \,sr} =({1}/{\rho_{\rm c}})({{\rm d}\rho_{\rm GW}}/{{\rm d}\ln f}) \sim F_{\rm sr}{\rho_{\rm BH}}/{\rho_{\rm c}}$, with $\rho_{\rm GW}$ and $\rho_{\rm BH}$ the GW and BH-mimicker energy densities. In the mass range $10^4{{\rm \mbox{--}}}10^7 M_\odot$ relevant for LISA, $\rho_{\rm BH}\sim {{\mathcal O}}( 10^4) M_{\odot}/{\rm Mpc}^3$, which gives $\Omega^{\rm LISA}_{\rm GW, \, sr} \sim 10^{-9}$. To estimate the background in the LIGO band, note that the background from ordinary BH binaries is $\newcommand{\e}{{\rm e}} \Omega_{\rm GW, \, bin}\sim \eta_{\rm GW} F_{\rm m} \rho_{\rm BH}/\rho_{\rm c}$, with $\newcommand{\e}{{\rm e}} \eta_{\rm GW}\sim {{\mathcal O}}(1\%)$ the GW emission efficiency for BH binaries [62], and $F_{\rm m} \sim {{\mathcal O}}(1\%)$ [55] the fraction of stellar-mass BHs in merging binaries. This gives $\newcommand{\e}{{\rm e}} \Omega_{\rm GW, \, sr}/\Omega_{\rm GW, \, bin}\sim F_{\rm sr}/(\eta_{\rm GW} F_{\rm m})\sim 10^2$, and because the LIGO O1 results imply $\Omega_{\rm GW, \, bin}\lesssim 10^{-9}{{\rm \mbox{--}}}10^{-8}$ [52, 63], we obtain $\Omega^{\rm LIGO}_{\rm GW, \, sr}\lesssim 10^{-7}{{\rm \mbox{--}}}10^{-6}$.

As shown in figure 3, LIGO/Virgo and LISA are also able to place model-independent constraints on the stochastic signal from exotic compact objects. At design sensitivity, LIGO/Virgo can detect or rule out any model with $t_0<10^{14} t_H$, whereas LISA can go as far as $t_0<10^{12}t_H$. In other words, objects across which light takes 10¹⁴ or less dynamical timescales to travel are ruled out. Finally, we note that although LIGO/Virgo rule out a wider range of the parameter space compared to LISA, it is still interesting to consider the constraints from both detectors since they probe different BH populations.

Zoom In Zoom Out Reset image size

Our results suggest that the current upper limits on the stochastic background from LIGO O1 already rule out the simplest models of BH mimickers at the Planck scale, setting the strongest constraints to date on exotic alternatives to BHs. The most relevant parameter for our analysis is the light travel time within the object, t₀. LIGO/Virgo (LISA) can potentially rule out models where $t_0<10^{14}t_H$ ($t_0<10^{12}t_H$). These results are not significantly dependent on the astrophysical uncertainties of the extragalactic BH distributions, and even if exotic compact objects are produced all at low spin, the order of magnitude of our constraints is unaffected.

We conclude that all present and future models of exotic compact objects should either conform with $t_0\gg 10^{14} t_H$ or represent at most a fraction X of the compact-object population, the remaining being BHs. Since equation (8) scales linearly with X, figure 2 implies that the O1 upper limits impose $X<50\%$ even if all BH-like objects are formed with low spins. At design sensitivity, these constraints could improve to $X<1\%$.

It might be possible to evade these constraints, by incorporating some (exotic and still unclear) mechanism quenching the ergoregion instability, e.g. absorption rates several orders of magnitude larger than those of neutron stars [35]. While such quenching mechanism might result in thermal or quasi-thermal electromagnetic radiation (which can be constrained by electromagnetic observations of BH candidates [64, 65]), quantum-dressed BH mimickers might evade such constraints by trapping thermal energy in their interiors for very long timescales [8, 15].

Our results also imply that in the simplest models of non-dissipative exotic ultracompact objects, GW echoes [6, 8, 25] can appear in the post-merger phase only after a delay time $\tau_{\rm echo}\sim t_0\gg 10^{14} t_H\approx 10^{10}[M/(20M_\odot)]{\rm s}$, which is orders of magnitude longer than what was claimed to be present in LIGO/Virgo data [66–68] (see also [69]). The latter claims would not be in tension with our bounds on the GW stochastic background only if one postulates exotic objects that are dissipative enough to absorb at least ${{\mathcal O}}(0.1)\%$ of gravitational radiation (which is several orders of magnitude more than what typically achievable with nuclear matter). In this case, the stochastic background from a population of 'echoing' merger remnants might still be detectable [70].

Finally, our constraints are stronger than those one might infer from BH spin measurements in x-ray binaries (see e.g. [71]). Indeed, those observations can only rule out instability timescales shorter than the BH age (in the more likely case in which the spin is natal [72], e.g. from supernova explosions) or the Salpeter timescale $t_{\rm S}\sim 10^7$ yr (if the BH spin is produced by accretion), otherwise the BH would have no time to spin down under the effect of the instability. However, the BH age is hard to reconstruct from the observed source properties, and in the few high-spin x-ray binaries where it was obtained via population synthesis modelling, it is of a few Myr [73, 74]. Therefore, obtaining exact bounds on the instability is tricky in the case of natal BH spins. Even if BH spins in x-ray binaries were accretion-produced, which is unlikely [72, 73], the obtained bound would simply be $t_{\rm inst}\gtrsim t_{\rm S}$, which would still be marginally weaker than the constraints presented in this paper (see figure 4).

Zoom In Zoom Out Reset image size

Let us also add that spin measurements in x-ray binaries are likely affected by unknown systematics (in several cases different techniques yield different results, see table 1 in [71]). More importantly, the very existence of the ergoregion instability in BHs surrounded by gas has never been investigated in detail, and the backreaction of the disk mass and angular momentum on the geometry, as well as the viscosity of the gas, may change the character and timescale of the instability.

Finally, let us comment on the relation between this work and [61], which computes the stochastic background from a hypothetical spin-down mechanism of BH merger remnants (in the LIGO/Virgo band only). Besides performing more realistic calculation of the spectrum ${\rm d}E/{\rm d}f$ due to the ergoregion instability, our work crucially differs from [61] as we compute the background from all BH mimickers in the Universe (including isolated ones), unlike [61] which only accounted for objects resulting from binary mergers. Since isolated compact objects are expected to be $1/F_{\rm m}\sim 100$ times more numerous than merging binaries [55], our background level is about 100 times larger, which allows constraining BH alternatives with current LIGO/Virgo data.

We thank Andrew Matas for a careful reading of the manuscript. We thank Jing Ren and Elisa Maggio for useful discussions. VC 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. PP acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement no. DarkGRA–757480. 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. The authors would like to acknowledge networking support by the COST Action CA16104.