All Dark Matter as Primordial Black Holes: a Gravitational-Wave Test

Introduction

The idea that black holes may have formed in the early Universe instead of as the final state of stellar evolution, goes back to the late sixties and early seventies [1, 2, 3, 4], and the suggestion that such primordial black holes (PBHs) could account for the dark matter (DM) without invoking any new elementary particle followed almost immediately [5, 6] (for recent reviews on the topic see Refs. [7, 8, 9, 10]).

The formation channels that have been proposed are numerous and predict very different mass spectra. The most studied one is the direct collapse of a large curvature perturbation generated during inflation once the corresponding scale re-enters the cosmological horizon, which, to allow a copious production of PBHs, requires an enhancement of the power spectrum with respect to the amplitude measured on cosmic microwave background scales. In practice, this can be realised either through a transient violation of the slow roll dynamics as the inflaton field evolves through its potential [11] or by the presence of a spectator field [12, 13, 14, 15, 16, 17, 18, 19, 20, 21].

Beyond the collapse of adiabatic perturbations, PBHs may be produced at a first-order phase transitions [22, 23, 24, 25, 26], from the collapse of cosmic string loops [27, 28, 29, 30] or of domain wall networks [31, 32, 33], or during an early matter-dominated or slow-reheating era, where the absence of pressure support makes the collapse far easier and the resulting holes are born with a sizeable spin [34, 35, 36, 37, 38, 39, 40].

As of today, independently of the formation channel, the PBHs are allowed to contribute to the totality of dark matter only in the so-called asteroidal mass window

Indeed, below roughly $10^{-16} M_\odot$ the Hawking evaporation time drops below the age of the Universe, so that the PBHs have either disappeared already or are still radiating enough to leave an imprint on the Galactic and extragalactic photon backgrounds [41]. Above roughly $10^{-11} M_\odot$, instead, the Einstein radius of the lens becomes large compared with the projected radius of the source star, finite-source effects no longer wash out the magnification, and the absence of microlensing events towards the Magellanic Clouds and M31 closes the window from the opposite side [42, 43].

However, the production of Gravitational Waves (GWs) could be a complementary tool to explore this mass range. In particular if the DM is composed solely of PBHs, then they would accumulate in the center of the galaxies, around astrophysical supermassive black holes (SMBHs) [44, 45, 46, 47] and fall inside it, in a coalescence process [48, 49]. This scenario is interesting since the relevant orbital frequencies are of the order of the millihertz, probed by LISA [50], and the mechanism that brings PBHs into their immediate neighbourhood is already part of the standard astrophysical dynamics.

A SMBH that grows adiabatically drags the surrounding DM inwards and leaves behind a steep density spike, whose slope and inner boundary have been computed in full general relativity for both the Schwarzschild and the Kerr geometry [46, 51, 52]. If the DM is made of PBHs, the spike is a granular population of compact objects on bound orbits, and those that end up close enough to the SMBH lose energy to GWs, spiralling inwards until they cross the innermost stable orbit and are swallowed [53, 54]. Each of these systems is an extreme mass-ratio inspiral with $q = m/M \sim10^{-22}$–$10^{-17}$, radiating in the millihertz band because the orbital frequency is fixed by $M$ alone. However, the single signals are far too weak to be resolved individually under any assumption, so that the real observable is the incoherent superposition of signals emitted by the entire population over all the galaxies in the Universe, in the same sense in which the stellar-mass objects inhabiting the same nuclei produce an astrophysical EMRI background [55, 56, 57].

The aim of this work is to compute the GW background generated by DM PBHs inspiralling into SMBHs and to determine for which SMBH masses and spins, PBH masses, and spike slopes the signal can be detected by LISA [50, 58] and DECIGO [59, 60].

This manuscript is organised as follows. In § 2 we define the DM spike and solve for its evolution under radiation reaction. In § 3 we compute the GW background for generic inclined orbits. In § 4 we focus on the equatorial case, where everything can be done analytically. In § 5 we discuss the detectability with LISA, and we conclude in § 6. Some technicalities are present in three appendices.

The DM spike and its evolution

Initial condition

We assume a spherically symmetric spike composed of PBHs with a power-law profile,

with $\rho_0$ the density at the outer radius $r_0$ and the slope in the range $1 \le\gamma< 3$, depending on the astrophysical processes that shape the spike [46, 51, 52, 54, 61, 62].

The lower limit $\gamma= 1$ corresponds to the inner cusp of a Navarro-Frenk-White halo in the absence of a central SMBH, while $\gamma\to3$ is the upper limit for a power-law spike of finite mass. Two-body relaxation of an equal-mass population of compact objects leads to the Bahcall-Wolf profile, with $\gamma= 7/4$ [63, 64]. The adiabatic growth of a SMBH steepens instead a pre-existing DM cusp into a spike with $2.25 \le\gamma\le2.5$, depending on the slope of the initial halo [46, 53, 54, 65].¹

We assume that the PBHs in the spike can be treated as collisionless particles moving along spherical geodesics around central SMBH, to simplify our calculations. Therefore we can set the inner boundary $r_{\rm in}$ as the innermost stable orbit. For equatorial orbits this is the usual innermost stable circular orbit (ISCO), while for inclined orbits it becomes the innermost stable spherical orbit (ISSO), which depends on the inclination angle $\iota$ and lies between the prograde and the retrograde equatorial values,

being $\chi$ the dimensionless spin. The orientation of the orbit is parametrised by $x \equiv\cos\iota$, where $x = 1$ and $x = -1$ correspond to equatorial prograde and retrograde orbits, respectively, and $x$ is a conserved quantity for the timelike geodesics, as explained in § 3.

For an isotropic distribution, the probability density in $x$ is uniform,

For a monochromatic mass function ² the number density is $n_{\rm PBH} = \rho_{\rm PBH}/m$, and the initial distribution of the PBHs in a single spike reads

where the factor $4\pi r^2$ follows from the spherical symmetry of the spike and the factor $P(x) = 1/2$ from the isotropy of the orbital orientations. We define the total mass of the spike as a fraction $\Delta$ of the SMBH mass, $\Delta\equiv M_{\rm sp}/M$, with

Throughout this work we fix the outer radius of the spike to $r_0 = 2GM/3\sigma_{\mathrm{DM}}^2$, using the relation between the SMBH mass $M$ and the fitted velocity dispersion $\sigma_{\mathrm{DM}}$ from Ref. [67]. For $\gamma< 3$ the term $r_{\mathrm{ISSO}}(\chi, x)$ is always negligible with respect to the outer edge $r_0$, so we can invert invert the relation above to get

The spike fraction $\Delta$ is not a free parameter of the model, since a spike massive enough to affect the stellar or gas kinematics would already have been absorbed into the dynamical mass inferred for the central black hole [68, 69]. Direct measurements are available only for a few galaxies, for example $\Delta\lesssim0.5$ for Centaurus A [70], $\Delta\lesssim0.08$ for Sgr A* [71] and $\Delta\lesssim0.1$ for M87 [72], so it can be taken of the order of 0.1. An alternative definition, discussed in Ref. [54], relates $\Delta$ to the slope of the spike through the influence radius of the SMBH [73]. The DM mass enclosed within a radius $R$ is

Defining the influence radius $r_h$ as the radius for which the DM mass is 2 times the central mass object [73] we have

and we set $r_0 = 0.2r_h$ to close the relation, as in [74]. Extending the profile to $R = r_h$ we obtain

that leads to

In the following, we use this relation to compute the amplitude of the signal.

Radiation reaction and the temporal evolution of the spike

For a quasi-circular orbit the radiation reaction integrated between $r$ and the ISSO yields the inspiral time [75, 76]

where in the last step we used $m\ll M$. Inverting this relation at fixed inclination gives the radius within which all the PBHs have reached the ISSO, and have therefore merged by the time $\tau(z)$,

where

is the cosmic time elapsed since the formation of the spike, which we set at the fiducial redshift $z_{\mathrm{form}}=5$ within the range of formation epochs considered in previous studies of cosmological DM spikes [77, 78]. In the range $0\leq z\leq z_{\mathrm{form}}$ we can write the elapsed time analytically as

where $H_0$ is the Hubble parameter and $\Omega_\Lambda$, $\Omega_{\mathrm{matter}}$ are the density fractions of the cosmological constant and matter respectively.

We now study how the spike evolves with time. A PBH released at radius $r_i$ is found at time $\tau$ at the radius $r$ fixed by $\tau=k(r_i^4-r^4)$, that is

Since the PBHs are neither created nor destroyed before reaching the inner boundary, their number in corresponding shells is conserved³,

Using the initial isotropic distribution

we obtain

on $r_{\mathrm{ISSO}} \le r_i(r,\tau) \le r_0(M)$. For $r^4 \gg\tau/k$ the spike can be considered frozen⁴, and Eq. (19) reduces to the initial profile (5).

The GW background from PBH spikes

Spherical orbits in Kerr geometry

We study now the orbits in a spinning BH background. The Kerr geometry can be written in Boyer-Lindquist coordinates as

where $M$ is the mass of the black hole, $aM$ is its angular momentum, and

A timelike geodesic is characterized by three constants of motion: the specific energy $E$, the specific axial angular momentum $L_z$ and the Carter constant $Q$ [79]. Therefore we define

In terms of the Mino time $\lambda$, the radial and polar motions decouple completely,

where $d/d\lambda=\Sigma\,d/d\tau$, with $\tau$ the proper time of the particle.

Orbits with zero eccentricity are spherical, which means that the radius is constant while the particle oscillates in the polar direction between two turning points and rotates in the azimuthal one. For a given radius and inclination these orbits are obtained by imposing $R(r)=R'(r)=0$, which fixes the energy $E$ and the axial angular momentum $L_z$ [80], while the additional condition $R''(r)=0$ defines the innermost stable spherical orbit $r_{\mathrm{ISSO}}(\chi,x)$. From the equations above we can extract the physical orbital frequencies $d\theta/dt=\Omega_\theta(r,x)$ and $d\phi/dt=\Omega_\phi(r,x)$ for every value of the inclination $x$ and of the spin $\chi$. We compute these frequencies and $r_{\mathrm{ISSO}}(\chi,x)$ with the KerrGeodesics package of the Black Hole Perturbation Toolkit [81, 82].

Energy flux to infinity

The GW flux emitted by a point particle moving on a bound Kerr geodesic is computed by solving the Teukolsky equation with a source for the perturbations of the Kerr spacetime [83, 84, 85]. Since the source motion is quasi-periodic, each mode can be labeled by the quantities $(s,\ell,\tilde{m},k,n)$, with $s=-2$ for gravitational perturbations, $\ell$ and $\tilde{m}$ the spheroidal-harmonic indices and $k$ and $n$ the polar and radial harmonic numbers (we denote the azimuthal number by $\tilde{m}$ to avoid confusion with the PBH mass $m$), and has frequency

We evaluate the modes with the TeukolskyPointParticleMode function of the Black Hole Perturbation Toolkit, keeping only

where $n=0$ because the orbits are spherical and $k=0$ is a conservative choice.⁵ The frequency of the selected modes is

where we restored the dimensions. Summing the two contributions $\tilde{m}=\pm\ell$, the luminosity at infinity of each multipole is

The inclination-averaged GW background

At fixed $x$ the relation (29) is monotonic and can be inverted numerically to give $r_{\ell}=r(f_{\ell},x)$, so that both the luminosity and the distribution of the PBHs can be written as functions of the frequency. The spectral luminosity of a single spike then follows by weighting the emission of each orbit with the number of PBHs that occupy it, summing over the multipoles and integrating over the orientations,

where the two Heaviside functions restrict the emission of each multipole to the frequencies spanned, at the given inclination, by the orbits between $r_{0}$ (which corresponds to $f_{\min}$) and the ISSO. The background observed today follows from [86], with the source frame frequency redshifted as $f_{s}=(1+z)f$,

where $\rho_c = 3H_0^2c^2/(8\pi G)$ is the critical energy density and $dn/d\ln M$ the comoving SMBH mass function. For the latter we adopt the optimistic scenario of Ref. [57], based on the occupation-corrected local black-hole mass function of Ref. [87], and we discuss the impact of this choice in § 5.

The equatorial limit

The equatorial prograde orbit, $x=1$, has the smallest ISCO radius and is therefore the most energetic orbital configuration, providing an upper estimate of the signal. The equatorial retrograde orbit, $x=-1$, has instead the largest ISCO radius and provides a natural lower bound.

Circular equatorial geodesics of the Kerr metric have angular frequency

so that the frequency of the dominant quadrupolar emission is

with the convention that $\chi>0$ and $\chi<0$ describe prograde and retrograde orbits around a black hole of spin $|\chi|$. Inverting the relation between frequency and radius gives

The ISCO radius is given by [88]

with

The emission extends from the frequency of the outermost orbits, at $r_0$, up to the ISCO frequency,

where

Every PBH is a source radiating at the monochromatic frequency (34), and the spectral luminosity of the spike around each SMBH is

where we retained only the equatorial prograde orbits, $x = 1$, which give the largest contribution, and where we remind that $f_s = (1+z)f$ and $P(x) = \frac{1}{2}$. We approximate the luminosity with the Newtonian quadrupole formula for circular orbits for simplicity, so we use

Including the evolution of the spike, described in Eq. (19)

we obtain a closed expression for the spectral luminosity

where we introduced $r_g = GM/c^2$ and restored $G$ and $c$. The background then follows from Eq. (32). This approximation allows us to study how the peak of the signal depends on the masses of the SMBH and the PBH and the parameters of the spike. We highlight the linear scaling with the PBH mass $m$ in the scenario when the contribution of $\tau(z)$ is negligible, i.e., in the frozen spike limit. Furthermore the factor $c/\pi r_g f_s-\chi$ shows that at higher spins the signal is boosted, given that the ISCO lies closer to the SMBH horizon.

Figure 2 compares the GW background obtained from the full numerical calculation of Eq. (32), which includes all the orbital inclinations and uses only the quadrupolar modes $\ell= 2$, $\tilde{m} = \pm2$, with the analytical approximation based only on equatorial prograde orbits. The analytical result captures the shape of the spectrum quite well but overestimates the peak by about one order of magnitude and shifts it slightly to higher frequencies. The difference is due to the Newtonian quadrupole luminosity, which is not reliable in the strong-field region, and we study this correction in App. § 9.

Prospects for future GW interferometers

In Figure 3 we show the GW background of Eq. (32) together with the power-law integrated sensitivities of LISA and DECIGO. All curves⁶ include the multipoles up to $\ell= 8$, with the modes $\tilde{m} = \pm\ell$. Our curves are obtained considering only the contribution coming from SMBH masses $5 \times10^5 M_\odot\le M \le10^7 M_\odot$ with the optimistic mass function of Ref. [57]. For these masses X-ray reflection measurements suggest very high spins [89], so we take $\chi= 0.95$ as fiducial value. Such a high spin is less justified outside this mass range, and restricting the integration to it gives a conservative estimate of the signal. The implications of the mass cut and an estimate of the GW signal in the entire mass range are discussed in App. § 11.

Dependence on the physical parameters

Spike slope. The top panel of Figure 3 shows how the signal depends on the slope of the spike. Our benchmark value is $\gamma= 2.5$, the steepest profile produced by the adiabatic growth of the SMBH [46]. We also show $\gamma= 2.75$ as a phenomenological example of a steeper profile [66]. The signal grows quickly with $\gamma$, because a steeper profile places more PBHs at small radii, where the orbital frequency and the emitted power are higher. Going from $\gamma= 2.5$ to $\gamma= 2.75$ raises the peak of the spectrum by an order of magnitude. As a consequence, for $\gamma= 2.75$ the central prediction enters the LISA sensitivity region, while for $\gamma= 2.5$ it remains slightly below the sensitivity curve. DECIGO, instead, can probe both configurations thanks to its lower sensitivity curve. The shaded bands show the uncertainty on the normalisation of the signal. Three effects can increase it:

• Close to a rapidly rotating black hole, the density of the spike computed in the Kerr geometry exceeds the Schwarzschild one by a factor of two to three [52].

• The outer radius of the spike can be identified with the Bondi radius $r_{\mathrm{B}} = 2GM/c_{\mathrm{s}}^2$, or with $r_0 = 2GM/v_{\mathrm{dm}}^2 = (5/9)r_{\mathrm{B}}$. The second choice is appropriate for a collisionless population, whose effective sound speed is $c_{\mathrm{s,eff}}^2 = (5/9)v_{\mathrm{dm}}^2$, with $v_{\mathrm{dm}} = 3\sigma_{\mathrm{DM}}$ [54, 74]. Since $\rho_0 r_0^\gamma\propto r_0^{\gamma-3}$, the smaller radius enhances the signal by $(5/9)^{\gamma-3}$.

• We neglect the modes with $\tilde{m} \ne\pm\ell$ or $k \ne0$. In App. § 10 we find that the polar modes $\ell= 2$, $\tilde{m} = 0$, $k = \pm2$ alone carry about one half of the power of the azimuthal modes $\ell= 2$, $\tilde{m} = \pm2$, $k = 0$.

The product of these three points, in the fiducial scenario $\gamma= 2.5$, is around $3 \times1.3 \times1.5 \simeq6$. In the opposite direction, the main uncertainty comes from the SMBH mass function. Replacing the mass function with the power-law fit [57, 93] $dn/d\ln M = 0.005 (M/4 \times10^{6} M_\odot)^{-0.3}\,\mathrm{Mpc}^{-3}$ lowers the signal in the LISA band by about a factor of four. The band extends from $\Omega_{\mathrm{GW}}/4$ to $6\Omega_{\mathrm{GW}}$.

PBH mass and SMBH spin. The bottom panel of Figure 3 fixes $\gamma= 2.5$ and compares $m = 10^{-11} M_\odot$ and $m = 10^{-12} M_\odot$, for our fiducial spin $\chi= 0.95$ and for $\chi= 0.8$, the weighted average of the spins reported in [89].

Close to the peak the growth with $m$ is slower, because heavier PBHs inspiral faster and have already emptied the innermost orbits, see Eq. (19). A larger spin brings the ISSO closer to the horizon, where the orbital frequency and the emitted power are higher. This raises the peak and extends the spectrum to higher frequencies. At low frequencies the signal is linear in $m$ and the emission comes from PBHs far from the SMBH. There the spin plays no role, and the curves for the two values of $\chi$ coincide. At the peak, lowering $\chi$ from 0.95 to 0.8 reduces the signal by approximately a factor of two. Within the range of parameters considered, the detectability of the signal depends more strongly on the spike slope rather than the SMBH spin. A detailed modeling of the spike slope is therefore essential to determine whether the GW background can be detected.

Signal-to-noise ratio analysis for LISA

To quantify the detectability with LISA we compute the signal-to-noise ratio (SNR) of the signal for a given slope and spin. The SNR is defined as

When foregrounds are absent the noise corresponds to the detector noise and thus we can write

and $S_h(f)$ is the LISA strain sensitivity [90, 91, 92]. Astrophysical foregrounds act as a further source of noise, reducing the sensitivity to stochastic signal. Different foregrounds are expected in the LISA band; in our estimates we include an extragalactic foreground from black-hole and neutron-star binaries [94, 95], a Galactic foreground from double white-dwarf binaries [96, 97, 98], and an extragalactic foreground from double white-dwarf binaries [99]. To account for the presence of foregrounds we include their spectra to the LISA noise in $\Omega_{\mathrm{noise}}$ in Eq. (44), following the analytical shapes provided by [100].

Figure 4 shows the SNR for $T_{\mathrm{obs}} = 4$ yr as a function of the slope $\gamma$ and of the PBH mass, up to the upper edge of the asteroidal window, $m = 10^{-11} M_\odot$. Since the signal is almost linear in $m$ where LISA is most sensitive, the SNR grows by about a factor of ten for each decade in $m$. The foregrounds lower the SNR by about 25%, so the contours with and without foregrounds almost overlap. For $\gamma\lesssim2.4$ the SNR remains below unity for every PBH mass in the window. Only sufficiently steep spikes place enough PBHs close to the central BH to produce a large SNR. There the orbital frequency is higher and the emitted power grows rapidly with decreasing radius, $L_{\mathrm{GW}} \propto r^{-5}$, so these PBHs dominate the signal. For the slope $\gamma= 2.5$ an SNR above one requires $m \gtrsim2 \times10^{-12} M_\odot$, and at $m = 10^{-11} M_\odot$ the SNR is about five. An SNR of ten requires $\gamma\gtrsim2.55$ for $m = 10^{-11} M_\odot$ and $\gamma\simeq2.7$ for $m = 10^{-12} M_\odot$, and it is not reached for $m \lesssim5 \times10^{-13} M_\odot$ within the slope range we consider.

Conclusions

In this work we have computed the stochastic GW background generated by asteroidal-mass PBHs that coalesce into the SMBH at the center of galaxies, after accumulating in the surrounding DM spikes. Although these systems have extremely small mass ratios and cannot be resolved individually, their combined emission over all galaxies produces a background in the millihertz band.

We have computed the GW background numerically for generic inclined spherical orbits in the Kerr geometry. We have also derived an analytical approximation based on equatorial prograde orbits, which have the smallest ISCO and give the largest contribution among all inclinations. This approximation reproduces the frequency dependence and peak amplitude of the numerical result.

The amplitude depends strongly on the spike slope. For our fiducial predictions, detection with LISA requires slopes above $\gamma= 2.5$, the largest value produced by adiabatic SMBH growth. Even a modest increase above this value can raise the GW background enough to make it detectable. Nevertheless, these predictions remain subject to the theoretical uncertainties of the model. A more complete calculation, including all the emission modes, a proper integration over the full SMBH mass range (that should involve a mass distribution depending on spin and redshift) and a spike profile computed in the Kerr geometry, could bring the signal within reach of LISA also for the value $\gamma= 2.5$. DECIGO can instead detect the signal for all the parameter configurations studied. In this case, a proper understanding of the spin distribution of low mass SMBHs would be needed, since it would affect the shape of the final spectrum.

This GW background provides a way to probe asteroidal-mass PBH dark matter independently of its formation mechanism and its detectability depends on the properties of the DM spikes around SMBHs, with the spike slope playing the main role in the scenarios considered.

Acknowledgments

A.R. and A.I. acknowledge support from the Swiss National Science Foundation (project number CRSII5_213497). A.J.I.’s research is funded by Tamkeen under the NYU Abu Dhabi Research Institute grant ADHPG-AD457. G.P. is supported by the Estonian Research Council grants PSG869, KOHTO34, TARISTU24-TK3, TARISTU24-TK10, and the Center of Excellence program TK202. D.P. is supported by the Swiss National Science Foundation under grants no. P500-2_239221. The Center of Gravity is a Center of Excellence funded by the Danish National Research Foundation under grant no. DNRF184.

Appendices

Relativistic correction to the flux

In this appendix we quantify the breakdown of the Newtonian approximation in the strong-field regime and use the relativistic correction to assess the accuracy of the analytical prediction for the GW background. The left panel of Figure 5 shows the ratio between the Newtonian quadrupole luminosity and the numerical Teukolsky flux for an equatorial prograde orbit, $x = 1$, considering the dominant modes $\ell= 2$, $\tilde{m} = \pm2$. At large radii the ratio approaches unity, as expected in the weak-field limit, while towards the ISCO the relativistic effects become important and the Newtonian luminosity exceeds the Teukolsky flux by about a factor of eight.

In the right panel we compare the GW background obtained from the full numerical calculation of Eq. (32), including all the orbital inclinations and using only the quadrupolar contributions $\ell= 2$, $\tilde{m} = \pm2$, with the analytical approximation based only on equatorial prograde orbits, where now, contrary to Figure 2, the analytical luminosity is corrected with the radius-dependent ratio shown in the left panel. The corrected analytical result reproduces in a good approximation the amplitude at the peak of the numerical spectrum, within a factor of about three. This result shows that the prograde configuration captures the dominant part of the emission once the luminosity is computed in the strong-field regime. At the same time, the comparison shows that a relativistic flux is needed for reliable predictions close to the ISCO, where the Newtonian approximation breaks down. The remaining discrepancy between the numerical and the analytical curve is expected, since the prograde approximation gives an upper estimate of the emission, whereas the numerical result averages over all orbital inclinations. We notice also that at low frequencies, the two curves coincide as the Newtonian limit is completely recovered.

Estimate of the polar-mode contribution

In the numerical calculation of the main text we retain the dominant azimuthal modes $\widetilde{m} = \pm\ell$ with $k = 0$. To estimate the contribution neglected by this truncation we compute separately the leading polar modes

which are relevant especially for polar orbits, $x \simeq0$. For spherical orbits $\Omega_r = 0$, and for $\widetilde{m} = 0$ the frequency of the mode is entirely determined by the polar motion. Restoring the dimensions, it reads

The two modes $k = \pm2$ have the same physical frequency once the absolute value is taken, and their fluxes add up to the total luminosity at infinity

At fixed inclination Eq. (B.2) can be inverted numerically to obtain $r_k = r(f_k,x)$, and the inclination-integrated spectral luminosity of the polar modes is

We compare this result with the inclination-integrated luminosity of the dominant quadrupolar azimuthal modes, $(\ell,\widetilde{m},k)=(2,\pm2,0)$, and at each frequency in their common domain we find

The polar modes alone contribute about one half of the dominant azimuthal component, which shows that the modes omitted from the calculation can increase the overall normalisation of the GW background and motivates their inclusion in the upper edge of the uncertainty band. We therefore conclude that the dominant energetic contribution comes from the quadrupolar azimuthal modes $\ell= 2$, $\tilde{m} = \pm2$ for equatorial prograde orbits, $x = 1$.

The spin and the full mass range

In the main text we integrated Eq. (32) only for SMBH masses in the range $5 \times10^{5}\,M_\odot\le M \le10^{7}\,M_\odot$. For these masses X-ray reflection measurements indicate high spins [89] that we have fixed correspondingly to $\chi= 0.95$. Outside this range the spin is poorly known, even though there is a tendency towards lower spins for lower masses. Since the spin fixes the position of the ISSO and is crucial to evaluate the resulting GW spectrum, in this appendix we estimate the GW signal in the entire mass range.

To do so, we extend the integration of Eq. (32) to the full range $10^{4}\,M_\odot\le M \le10^{9}\,M_\odot$, considering $\chi= 0.95$ and $\chi= 0.5$. For both missing ranges, i.e. for SMBHs of masses $10^{4}\,M_\odot\le M \le5 \times10^{5}\,M_\odot$ and of masses $10^{7}\,M_\odot\le M \le10^{9}\,M_\odot$, we estimate the GW spectrum produced for two different values of the spin, based on [89]: a maximal one $\chi= 0.95$ and a lower one $\chi= 0.5$. The spectra evaluated at these two values determine the band in which the effective spectrum is expected. In this way, not only we estimate the impact of those masses on the final spectrum, but also the corresponding frequency region which would be more affected.

The results are shown in Figure 6. In red we show our benchmark signal, obtained for $5 \times10^{5}\,M_\odot\le M \le10^{7}\,M_\odot$ and $\chi= 0.95$. The band corresponding to higher masses is reported in yellow. These masses would mostly contribute to the IR tail, with a quite negligible contribution in the range of interest. The band for lower masses, instead, is reported in blue. Those would dominate the signal at higher frequencies and have a major effect observable by DECIGO. We checked that, even in the most optimistic case, these contributions would leave the portion of the spectrum that falls within the LISA sensitivity band practically unchanged. Neglecting them therefore does not affect the analysis we report in § 5.2. The situation is different for DECIGO, whose sensitivity covers the frequencies where the light SMBHs dominate. There the band gives the expected range of the signal, while a model of the spin distribution of light SMBHs would be needed to predict its detailed shape. For this reason in § 5.2 we limited our analysis to the LISA case. For reference, we also report in black the total signal that would be observed if the spin of all the SMBH in the whole mass range would follow a monochromatic distribution with $\chi= 0.95$.