Open Access
Issue
A&A
Volume 676, August 2023
Article Number A38
Number of page(s) 13
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202346268
Published online 02 August 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The discovery (LIGO Scientific & Virgo Collaborations 2016) of gravitational waves (GWs) by the LIGO and Virgo experiments, in the 𝒪(100) Hz range of frequencies generated by mergers of black holes (BHs) with masses of 𝒪(10 − 100) M (LIGO Scientific, Virgo & KAGRA Collaborations 2021), has opened a new window onto both astrophysics and cosmology. Supermassive black holes (SMBHs) with masses of 𝒪(106 − 1010) M are known to be present in galactic nuclei, and the immediate surroundings of two of them have recently been imaged by the Event Horizon Telescope (Event Horizon Telescope Collaboration 2019, 2022). However, information on intermediate-mass black holes (IMBHs) with masses in the range 𝒪(103 − 105) M (Greene et al. 2020) is less complete. Observations of IMBHs and their mergers would cast light on the uncharted mechanisms that must be presented for SMBHs to form.

A prerequisite for BH mergers is the formation of tightly bound binaries that can radiate GWs efficiently. While there are several well-understood paths for forming compact stellar-mass BHs binaries (e.g., common envelope evolution or dynamical capture in dense stellar environments), the final stage of SMBH binary formation in galaxy mergers is still not fully understood and is commonly known as “the final parsec problem” (Begelman et al. 1980). However, observations of quasar optical variability provide indirect evidence for the existence of some tight binary SMBHs in the regime where the emission of GWs must already have a noticeable effects on their orbital evolution (Rieger & Mannheim 2000; De Paolis et al. 2002, 2004; Valtonen 2008; Boroson & Lauer 2009; Iguchi et al. 2010; Graham et al. 2015; O’Neill et al. 2022; Kovačević et al. 2022)1. The GW data will be essential for a better understanding how SMBHs overcome the “final parsec” obstacle, and how the assembly of the SMBH population proceeds in general.

Pulsar timing arrays (PTAs) are potentially sensitive to GWs in the nHz range, while NANOGrav and other PTAs have recently reported evidence for a common-spectrum stochastic process (Arzoumanian et al. 2020; Goncharov et al. 2021; Chen et al. 2021; Antoniadis et al. 2022). Their signals have power spectra that are consistent with predictions based on inspiralling binary SMBH models (Phinney 2001)2, but they have not (yet) detected the Hellings-Downs quadrupolar signature that is characteristic of GWs (Hellings & Downs 1983)3.

The BH merger interpretation of the PTA measurements does not provide direct information on the masses of the infalling BHs, which requires the modelling and measurements of this stochastic GW background in different frequency ranges (Sesana et al. 2008). The purpose of this paper is to illustrate how future measurements of the GW spectrum due to unresolved infall sources and individual binary mergers at frequencies between the PTAs and LIGO/Virgo will be able to extend the PTA measurements to lower black-hole mass ranges, probing models for the assembly of SMBHs, and perhaps obtaining direct evidence for the mergers of IMBHs. We take as the amplitude of the common-spectrum process measured by the PTAs A = 2 . 8 0.8 + 1.2 × 10 15 $ A = 2.8^{+1.2}_{-0.8} \times 10^{-15} $ at a reference frequency of 1/yr, as found by the International Pulsar Timing Array (IPTA; Antoniadis et al. 2022) assuming the spectral index of α = −2/3 expected from inspiralling SMBHs, which is consistent with other measurements.

We use examples of next-generation GW detectors the planned LISA laser interferometer (Berti et al. 2006), whose peak sensitivity is at frequencies 𝒪(10−4 − 10−2) Hz, the proposed ET laser interferometer (Sathyaprakash et al. 2012), whose peak sensitivity is at frequencies 𝒪(1 − 102) Hz and the projected AION-km and AEDGE atom interferometers (Badurina et al. 2020; El-Neaj et al. 2020), whose peak sensitivities would be in the range of 𝒪(10−2 − 1) Hz. We show that, while none of these detectors would observe a signal from the binary black holes weighing > 𝒪(109) M that would probably be responsible for most of the PTA signal, LISA would observe a merger signal from binary black holes weighing 𝒪(103 − 109) M that might contribute part of the PTA signal, while AEDGE would observe a merger signal if there is a population of binary black holes weighing 𝒪(103 − 106) M with a formation history similar to the heavier BHs responsible for the PTA signal.

Throughout this paper, we use natural units with c = 1, GN = 1.

2. BH merger rate

As a first step towards estimating the BH merger rate, we use the Extended Press-Schechter (EPS) formalism (Press & Schechter 1974; Bond et al. 1991) to calculate the galactic halo mass function,

d n ( M , t ) d ln M = ρ 0 M 2 π d ln σ d ln M δ c ( z ) σ ( M ) e δ c ( z ) 2 2 σ 2 ( M ) , $$ \begin{aligned} \frac{\mathrm{d}n(M,t)}{\mathrm{d}\ln M} = \frac{\rho _0}{M} \sqrt{\frac{2}{\pi }} \frac{\mathrm{d}\ln \sigma }{\mathrm{d}\ln M} \frac{\delta _c(z)}{\sigma (M)}e^{-\frac{\delta _c(z)^2}{2\sigma ^2(M)}} , \end{aligned} $$(1)

where ρ0 is the background matter density today, σ2(M) is the variance of the matter fluctuations and δc(z) is the critical overdensity for collapse. The latter is given by δc(z)≈1.686/D(z) where D(z) is the linear growth function (Dodelson 2003). We calculate σ2(M) from the cold dark matter power spectrum with the Planck 2018 cosmological parameters (Planck Collaboration VI 2020) and the transfer function derived in Eisenstein & Hu (1998). The EPS formalism also gives an estimate for the probability per unit time for a halo of mass M1 to merge with another one of mass, M2, at some redshift, z, and become a halo of mass Mf = M1 + M2 (Lacey & Cole 1993):

d p ( M 1 , M 2 , t ) d t d M 2 = 1 M f 2 π | δ ˙ c δ c | d ln σ d ln M f [ 1 σ 2 ( M f ) σ 2 ( M 1 ) ] 3 2 × δ c σ ( M f ) exp [ δ c 2 2 ( 1 σ 2 ( M f ) 1 σ 2 ( M 1 ) ) ] . $$ \begin{aligned} \frac{\mathrm{d}p(M_1,M_2,t)}{\mathrm{d}t \mathrm{d}M_2}&= \frac{1}{M_f}\sqrt{\frac{2}{\pi }}\left|\frac{\dot{\delta }_c}{\delta _c}\right| \frac{\mathrm{d}\ln {\sigma }}{\mathrm{d}\ln {M_f}} \left[1-\frac{\sigma ^2(M_f)}{\sigma ^2(M_1)}\right]^{-\frac{3}{2}} \nonumber \\&\quad \times \frac{\delta _c}{\sigma (M_f)} \exp \!\left[-\frac{\delta _c^2}{2}\left(\frac{1}{\sigma ^2(M_f)}-\frac{1}{\sigma ^2(M_1)}\right)\right]. \end{aligned} $$(2)

Defining the halo merger rate kernel:

Q ( M 1 , M 2 , t ) d p ( M 1 , M 2 , t ) d t d M 2 [ d n ( M 2 , t ) d M 2 ] 1 , $$ \begin{aligned} Q(M_1,M_2,t) \equiv \frac{\mathrm{d}p(M_1,M_2,t)}{\mathrm{d}t \mathrm{d}M_2} \left[\frac{\mathrm{d}n(M_2,t)}{\mathrm{d}M_2} \right]^{-1} , \end{aligned} $$(3)

the halo merger rate is

d R h d M 1 d M 2 = d n ( M 1 , t ) d M 1 d n ( M 2 , t ) d M 2 Q ( M 1 , M 2 , t ) . $$ \begin{aligned} \frac{\mathrm{d}R_h}{\mathrm{d}M_1 \mathrm{d}M_2} = \frac{\mathrm{d}n(M_1,t)}{\mathrm{d}M_1} \frac{\mathrm{d}n(M_2,t)}{\mathrm{d}M_2} Q(M_1,M_2,t) . \end{aligned} $$(4)

As discussed by Benson et al. (2005), the merger rate kernel Q(M1, M2, t) in the EPS formalism is not symmetric in the exchange of M1 and M2. However, it was shown by Erickcek et al. (2006) that taking M1 < M2 in the merger kernel aptly approximates the result of the Benson-Kamionkowski-Hassani (BKH) merger theory (Benson et al. 2005) where the merger rate kernel is computed from the Smoluchowski coagulation equation and preserves the EPS halo mass function. We will, therefore, require that the first argument of Q is smaller than the second.

In order to estimate the merger rate of central BHs, we must estimate the probability pocc(m|M, z) that a BH of mass m occupies a halo of mass M, as well as the probability pmerg(m1, m2) that the galactic merger leads to a merger of their central BHs. The resulting BH merger rate can be expressed as

d R BH d m 1 d m 2 = d M 1 d M 2 p merg ( m 1 , m 2 ) × p occ ( m 1 | M 1 , z ) p occ ( m 2 | M 2 , z ) d R h d M 1 d M 2 . $$ \begin{aligned}&\frac{\mathrm{d}R_{\rm BH}}{\mathrm{d}m_1 \mathrm{d}m_2} = \int \mathrm{d}M_1 \mathrm{d}M_2 \, p_{\rm merg}(m_1,m_2) \nonumber \\&\qquad \qquad \,\,\, \times p_{\rm occ}(m_1|M_1,z) p_{\rm occ}(m_2|M_2,z) \frac{\mathrm{d}R_h}{\mathrm{d}M_1 \mathrm{d}M_2} . \end{aligned} $$(5)

We assume the following simple redshift-dependent relation between the halo mass and the BH mass (Barkana & Loeb 2001; Wyithe & Loeb 2003):

M v 10 12 M = 10.5 [ Ω M ( 0 ) Ω M ( z ) Δ c ( z ) 18 π 2 ] 1 2 ( 1 + z ) 3 2 [ m BH 10 8 M ] 3 5 , $$ \begin{aligned} \frac{M_v}{10^{12}M_\odot } \!=\! 10.5 \left[ \frac{\Omega _M(0)}{\Omega _M(z)} \frac{\Delta _c(z)}{18\pi ^2} \right]^{\!-\frac{1}{2}} \!\!(1+z)^{-\frac{3}{2}} \!\left[\frac{m_{\rm BH}}{10^8 M_\odot }\right]^{\!\frac{3}{5}}, \end{aligned} $$(6)

where Δc(z) = 18π2 + 82[ΩM(z)−1] − 39[ΩM(z)−1]2 is the critical overdensity at virialization, which corresponds to pocc(m|M) = pocc(m)δ(m − mBH(M)). It is generally expected that SMBHs inhabit most large galaxies, and X-ray observations (Miller et al. 2015) constrain the SMBH occupation fraction pocc(m) to be > 20% for early galaxies with lower stellar masses 107 < M*/M < 1010. It has been estimated that an initially small occupation fraction may grow at low redshifts (see, e.g., Lippai et al. 2009), which we neglect in this analysis. As we show later, the PTA GW signal is dominated by events in a range 1 ≲ z ≲ 3. The estimated z dependence of pocc would shift this range towards lower redshifts when compared to a constant pocc, increasing the observability of BH mergers.

In view of the paucity of information about the merger probability pmerg(m1, m2), we model pmerg by a constant for the sake of simplicity. After the galactic merger, various dynamical mechanisms must decrease the size of the SMBH binary below sub-parsec scales in order for a GW emission-driven merger to take place (Begelman et al. 1980). The crossing of the final parsec is determined mostly by binary hardening via stellar loss-cone scattering (Merritt 2013), which can currently be considered the largest source of uncertainty in pmerg. The numerical study of Kelley et al. (2017a), based on a population of 106 − 1010M BH binaries in the Illustris simulation (Vogelsberger et al. 2014), found for a wide range of model parameters that the coalescing fraction is nearly independent of the total mass of the binary but decreases with the mass ratio. Neglecting the latter dependence should not significantly affect our results, since extreme mass ratio inspirals contribute subdominantly to the GW background. With these assumptions, we obtain:

d R BH d m 1 d m 2 p BH d M 1 d m 1 d M 2 d m 2 d R h d M 1 d M 2 , $$ \begin{aligned} \frac{\mathrm{d}R_{\rm BH}}{\mathrm{d}m_1 \mathrm{d}m_2} \approx p_{\rm BH} \frac{\mathrm{d}M_1}{\mathrm{d}m_1} \frac{\mathrm{d}M_2}{\mathrm{d}m_2} \frac{\mathrm{d}R_h}{\mathrm{d}M_1 \mathrm{d}M_2} , \end{aligned} $$(7)

where pBH ≡ pocc(m1)pocc(m2)pmerg denotes the probability that the halo merger leads to a merger of the central SMBHs.

We show in Fig. 1 the BH merger rates as functions of the redshift z, calculated using this approach in three ranges of the total merging BH mass M and normalized relative to the merger probability pBH discussed above. We see that mergers with total masses M > 109M occur typically at z = 𝒪(1), those with total masses M ∈ (106, 109) M occur typically at z = 𝒪(3), and those with total masses M ∈ (103, 106) M occur typically at z = 𝒪(5).

thumbnail Fig. 1.

BH merger rates calculated as functions of the redshift z in three ranges of the total merging BH mass M and normalized relative to the merger probability pBH discussed in the text.

In this paper, we use the IPTA measurement (Antoniadis et al. 2022) to normalize pBH for large BH masses, which circumvents the astrophysical uncertainties related to pocc and pmerg, and extrapolate it to smaller BH masses. We also comment on the potential mass dependence of pBH, which is relevant for the GW phenomenology when we extrapolate from the frequency band relevant for PTAs to those to be explored by next-generation GW detectors.

3. Analysis

3.1. The GW energy spectrum

The mean GW energy density spectrum from the SMBH binary population can be estimated as (Phinney 2001):

Ω GW ( f ) 1 ρ c d ρ GW d ln f = d λ 2 π 5 f 3 | h ( f ) | 2 ρ c , $$ \begin{aligned} \Omega _{\rm GW}(f) \equiv \frac{1}{\rho _{\rm c}}\frac{\mathrm{d}\rho _{\rm GW}}{\mathrm{d}\ln f} = \int \text{ d}\lambda \,\frac{2\pi }{5} \frac{f^3 |\tilde{h}(f)|^2}{\rho _c} , \end{aligned} $$(8)

where | h ( f ) | $ |\tilde h(f)| $ denotes the optimal amplitude of the Fourier transform of the GW strain, ρ c =3 H 0 2 /8π $ \rho_{\rm c} = 3H_0^2/8\pi $,

d λ = d m 1 d m 2 d z 1 + z d V c d z d R BH ( t ) d m 1 d m 2 , $$ \begin{aligned} \mathrm{d}\lambda = \mathrm{d}m_1 \mathrm{d}m_2 \,\frac{\mathrm{d}z}{1+z} \frac{\mathrm{d}V_c}{\mathrm{d}z} \frac{\mathrm{d}R_{\rm BH}(t)}{\mathrm{d}m_1 \mathrm{d}m_2} , \end{aligned} $$(9)

is the differential merger rate and Vc denotes the comoving volume available at redshift z. We compute | h ( f ) | $ |\tilde h(f)| $ using the inspiral-merger-ringdown template (Ajith et al. 2008):

| h ( f ) | = 5 24 M z 5 6 π 2 3 D L × { f 7 6 f < f merg f merg 1 2 f 2 3 f merg f < f ring f merg 1 2 f ring 2 3 σ 2 4 ( f f ring ) 2 + σ 2 f ring f < f cut , $$ \begin{aligned}&|\tilde{h}(f)| = \sqrt{\frac{5}{24}}\frac{\mathcal{M} _z^{\frac{5}{6}}}{\pi ^{\frac{2}{3}}D_L}\nonumber \\&\qquad \quad \,\,\,\times {\left\{ \begin{array}{ll} f^{-\frac{7}{6}}&f<f_{\rm merg} \\ f_{\rm merg}^{-\frac{1}{2}}f^{-\frac{2}{3}}&f_{\rm merg}\le f< f_{\rm ring} \\ f_{\rm merg}^{-\frac{1}{2}}f_{\rm ring}^{-\frac{2}{3}}\frac{\sigma ^2}{4(f-f_{\rm ring})^2+\sigma ^2}&f_{\rm ring}\le f < f_{\rm cut} , \end{array}\right.} \end{aligned} $$(10)

where the redshift-dependent chirp mass, ℳz, is given in terms of the total mass, M, and the symmetric mass ratio η ≡ m1m2/M2 of the binary by M z ( 1 + z ) M 0 = ( 1 + z ) M η 3 5 $ \mathcal{M}_z \equiv (1+z) \mathcal{M}_0 = (1+z) M \eta^{\frac35} $, and DL denotes the luminosity distance of the binary. This implies the canonical Ω GW f 2 3 $ \Omega_{\mathrm{GW}} \propto f^{\frac23} $ scaling during the inspiral phase (Phinney 2001). The frequencies fmerg, fring, fcut, and σ are parameterised as f j = η 5 3 ( a j η 2 + b j η + c j ) / ( π M z ) $ f_j = \eta^{\frac53}(a_j \eta^2 + b_j \eta + c_j)/(\pi \mathcal{M}_z) $ where aj, bj, and cj are coefficients whose fitted values are given in Table I of Ajith et al. (2008).

We comment briefly on potential uncertainties in the GW energy density spectrum. If orbital decay is driven by other processes in addition to GW emission, such as viscous drag, then the coalescence time is shortened and the total GW spectrum is suppressed by a factor of Ttot/TGW, where Ttot, TGW denote the characteristic hardening timescales (Kocsis & Sesana 2011; Kelley et al. 2017a). Here we omit this effect, assuming that the signal is dominated by binaries for which the hardening is driven mainly by GW emission. Additionally, we consider only circular binaries, which, in the inspiral phase, radiate monochromatically at twice the orbital frequency. We note, however, that the GW spectrum of eccentric binaries would contain higher harmonics, with the fundamental harmonic (at the orbital frequency) becoming dominant at eccentricities e > 0.4 (Taylor et al. 2016). In addition, the GW luminosity would be enhanced by a factor of ( 1 + 73 / 24 e 2 + 37 / 96 e 4 ) / ( 1 e 2 ) 7 2 $ (1 +73/24 e^2 + 37/96 e^4)/(1-e^2)^{\frac72} $ (Peters & Mathews 1963; Enoki & Nagashima 2007; Kelley et al. 2017b).

Figure 2 shows the GW energy density spectrum from BH binaries in different mass ranges. The black dashed curve corresponds to pBH = 1, namely, the assumption that each galactic merger would produce a BH merger with m2 > m1 > 103M. This naive assumption leads to a spectrum exceeding that observed in the PTA band.

thumbnail Fig. 2.

Mean GW energy density spectrum from massive BH mergers compared with the sensitivities of different experiments. The black dashed curve shows the case where all mergers of galaxies produce a GW signal for BHs heavier than 103M (pBH = 1). The colored bands show the spectra from SMBHs heavier than 109M (dark blue), from SMBHs in the range (106M, 109M); light blue), and from IMBHs with masses in the range 103 − 106M (green), assuming a universal efficiency factor p BH = 0 . 17 0.08 + 0.18 $ p_{\mathrm{BH}} = 0.17^{+0.18}_{-0.08} $. The shaded regions show the prospective sensitivities of IPTA30 (Chen et al. 2017; Kaiser & McWilliams 2021), LISA (LISA Collaboration 2017), AEDGE (El-Neaj et al. 2020; Badurina et al. 2021), AION-km (Badurina et al. 2020, 2021) and ET (Sathyaprakash et al. 2012).

We stress that the spectrum in Eq. (8), which is displayed in Fig. 2, does not necessarily correspond to a stochastic GW background as it does not distinguish between individual resolvable events and unresolvable events that would contribute to a GW background. The Ω GW f 2 3 $ \Omega_{\mathrm{GW}} \propto f^{\frac23} $ tail at low frequencies arises from nearly monochromatic signals generated by inspiralling BH binaries. Because of this, ΩGW will fluctuate depending on the specific realization of the binary population, so the spectrum in Fig. 2 should be interpreted as the mean ΩGW obtained by averaging over many potential realizations of the binary population. We return to the probability distribution of ΩGW in the next subsection and in Appendix A.

The grey band in Fig. 2 was obtained by choosing a fixed value of pBH so that the amplitude of the total GW spectrum from BH binaries matches the IPTA observations in the nHz range, which requires p BH = 0 . 17 0.08 + 0.18 $ p_{\mathrm{BH}} = 0.17^{+0.18}_{-0.08} $. The colored bands show the contributions to the total IPTA spectrum from different ranges of the BH masses. In comparing the dark blue, light blue, and green bands, we can see that the largest contribution to the GW spectrum in the PTA window comes from M > 109M SMBH binaries, consistent with earlier studies (Wyithe & Loeb 2003; Sesana et al. 2004, 2008, 2009; Enoki et al. 2004; Kelley et al. 2017a; Izquierdo-Villalba et al. 2021; Bécsy et al. 2022).

Some characteristics of these binaries are shown in Fig. 3. The left panel shows that the total masses are ≲1010, and the redshifts are typically 𝒪(1). The middle left panel shows that the typical masses increase with decreasing frequency and the middle right panel shows that the binary lifetimes are only weakly correlated with their emission frequencies, with lifetimes ≳104 yr being generally favored. The right panel shows the GW strain from these binaries, together with the upper limit on continuous GW sources from NANOGrav 12.5 years data (Arzoumanian et al. 2023) and the projected sensitivity of IPTA30 (Chen et al. 2017; Kaiser & McWilliams 2021) at a signal-to-noise ratio of SNR = 8, which is similar to that expected for SKA (Janssen et al. 2015).

thumbnail Fig. 3.

Distributions of the binary parameters for a sampling of the BH binary population with 0 < z < 3, 1 nHz < f < 30nHz and M 0 > 10 9 ( f / nHz ) 4 3 $ \mathcal{M}_0 > 10^9 (f/\mathrm{nHz})^{-\frac43} $. The symmetric mass ratio is integrated with the lower bound η > 0.001. In the right panel, the black curve shows the NANOGrav 95% CL constraint on continuous GW sources (Arzoumanian et al. 2023) and the projected sensitivity of IPTA30 at SNR = 8.

Returning to Fig. 2, we see that the GW energy spectrum in the LISA frequency range is dominated by lighter BH binaries in the mass range of M ∈ (106, 109) M, shown in light blue, and the contribution in the AEDGE band is dominated by IMBHs with masses in the range of M ∈ (103, 106) M, shown in green. The low-frequency GW emissions from the early stages of the inspirals of these binaries also contribute in the lower-frequency bands, but these contributions are subdominant compared to the signals from higher-mass SMBH mergers4.

We recall that the estimated spectrum in Fig. 2 assumes a constant pBH. Clearly, the relative contribution from binaries with lighter BH masses < 109M could be enhanced (reduced) by increasing (decreasing) pBH for these binaries. In order for mergers of BHs with masses in the range (106, 109) M to dominate the PTA signal pBH ≃ 1 would be required in this mass range, but even this maximal enhancement would be insufficient for mergers with M < 106M to contribute significantly to the PTA signal. That said, we emphasize that the GW spectrum in either the LISA or AEDGE frequency range could reach the black dashed line if pBH ≃ 1 in the relevant mass range, and could also be enhanced if Population III stars make a significant contribution to the BH spectrum. Conversely, the spectrum could be suppressed in these ranges if pBH is smaller than the value that fits the PTA data.

3.2. Distribution of sources

The number of BH binaries at redshift z emitting in a given frequency band can be estimated from the time the binary spends in that frequency band, which gives (Sesana et al. 2008)

d N d M z d η d z d ln f = d λ d M z d η d z | d τ d ln f | , $$ \begin{aligned} \frac{\mathrm{d}N}{\mathrm{d}\mathcal{M} _z \mathrm{d}\eta \mathrm{d}z \mathrm{d}\ln f} = \frac{\mathrm{d}\lambda }{\mathrm{d}\mathcal{M} _z \mathrm{d}\eta \mathrm{d}z} \left| \frac{\mathrm{d}\tau }{\mathrm{d}\ln f} \right| , \end{aligned} $$(11)

where, assuming circular orbits and orbital decay via GW emission (Peters 1964), the coalescence time of an inspiralling binary emitting GWs with frequency (1 + z)f is5

τ ( f ) = 5 256 π 8 3 M z 5 3 f 8 3 . $$ \begin{aligned} \tau (f) = \frac{5}{256 \pi ^{\frac{8}{3}} \mathcal{M} _z^{\frac{5}{3}}} f^{-\frac{8}{3}} . \end{aligned} $$(12)

We show in Fig. 4 the expected numbers of binaries with total masses of M ≥ Mmin and mass ratios q ≥ 0.01 that emit GWs in several frequency bands in the range 1 nHz < f < 30 nHz (corresponding to the sensitivity ranges of PTAs). As seen in Fig. 2, SMBH binaries generate the dominant contribution to the GW signal in this frequency range. The fractions of the total GW signal in each frequency band that are generated by binaries with larger masses are indicated by the vertical dashed lines. For example, we see that almost 105 (102) binaries with M > 3 × 109M generate 80% of the calculated GW signal in the 1 nHz < f < 3 nHz (10 nHz < f < 30 nHz) frequency band, whereas 50% of the calculated GW signal in the 3 nHz < f < 10 nHz band may be generated by just one source. These examples indicate that the expected signal, particularly in higher-frequency bins, may be comprised of a limited number of nearly monochromatic signals from heavy SMBH binaries ( > 109M), implying sizeable fluctuations around the smooth f 2 3 $ f^{\frac23} $ spectrum that would be obtained in the limit of a large population of heavy inspiralling binaries.

thumbnail Fig. 4.

Expected number of SMBH binaries heavier than Mmin and with symmetric mass ratio η > 0.01 emitting GWs in the indicated frequency range for pBH = 0.17. The vertical dashed lines indicate the fractions of the total GW signal in each frequency band generated by binaries with larger masses.

In order to study the prospective SMBH binary population in the PTA band via a Monte Carlo approach, we generated realizations of the population from the probability density function of binaries emitting at a given frequency. We divide the frequency range into bins (fj, fj + 1) populated with N(fj) binaries. A realization of the GW background from each bin can then be obtained as (see Eq. (A.6)):

Ω GW ( f j ) = 1 ln ( f j + 1 / f j ) k = 1 N ( f j ) Ω GW ( 1 ) ( θ k ) , $$ \begin{aligned} \Omega _{\rm GW}(f_j) = \frac{1}{\ln (f_{j+1}/f_j)}\sum _{k=1}^{N(f_j)} \Omega _{\rm GW}^{(1)}(\boldsymbol{\theta }_k) , \end{aligned} $$(13)

where the contribution from an individual binary emitting in this frequency band is given by

Ω GW ( 1 ) ( θ ) = | d τ d ln f | 1 2 π 5 f 3 | h ( f ) | 2 ρ c , $$ \begin{aligned} \Omega _{\rm GW}^{(1)}(\boldsymbol{\theta }) = \left| \frac{\mathrm{d}\tau }{\mathrm{d}\ln f} \right|^{-1} \frac{2\pi }{5} \frac{f^3 |\tilde{h}(f)|^2}{\rho _c} , \end{aligned} $$(14)

and θ ≡ {ℳz,,f} are parameters describing the binary. The number of binaries N(fj) in this frequency bin is drawn from a Poisson distribution with the expectation value N ¯ ( f j ) f j f j + 1 d N $ \bar N(f_j) \equiv \int_{f_{j}}^{f_{j+1}} \!{\text{ d}}N $ and the binary parameters are generated randomly according to the distribution

p ( M z , η , z , f ) d N d M z d η d z d ln f . $$ \begin{aligned} p(\mathcal{M} _z, \eta , z, f) \propto \frac{\mathrm{d}N}{\mathrm{d}\mathcal{M} _z \mathrm{d}\eta \mathrm{d}z \mathrm{d}\ln f} \, . \end{aligned} $$(15)

This is the expression that was used to calculate the distributions of binary parameters shown in Fig. 3 for a sampling with pBH = 0.17 of 4 × 105 SMBH binaries with redshifts 0 < z < 3, mass ratios η > 0.01, and GW frequencies 1 nHz < f < 30 nHz. We display only binaries with chirp mass M 0 > 10 9 ( f / nHz ) 4 3 $ \mathcal{M}_0 > 10^9 (f/\mathrm{nHz})^{-\frac43} $, so that N ¯ ( f j ) > O ( 10 4 ) $ \bar N(f_j) > \mathcal{O}(10^4) $ for each 1 nHz frequency bin.

In order to generate the statistical distributions of ΩGW, it is sufficient to consider a smaller set of parameters, and the distribution Eq. (11) for the full population can be reduced to a simpler distribution P(1), f) for individual sources Ω GW ( 1 ) $ \Omega_{\mathrm{GW}}^{(1)} $ emitting at a frequency f, which can be expressed via two one-parameter functions depending on the merger rate model, as discussed in Appendix A.

The upper left panel of Fig. 5 illustrates the frequency spectra found in 1600 Monte Carlo realizations of the SMBH binary population drawn from a statistical distribution similar to that shown in Fig. 3. The black solid curve shows the mean of the spectra, the grey bands show the 1σ and 2σ CL regions of the GW spectra, and the red line is the spectrum found in one representative realization. We note the importance of fluctuations, particularly at higher frequencies where fewer SMBH binaries contribute, with individual binaries becoming distinguishable at frequencies f ≳ 10 nHz. We find that the median spectrum, shown by the dashed line, has a somewhat lower slope than the analytic result shown as the green line. Although the numerically obtained mean spectrum (black), lies below the analytic expectation, it will approach it slowly if the number of Monte Carlo realizations is increased. The tendency of the mean spectrum to lie along the upper side of the 1σ CL range is due to the fact that the spectrum has a long tail at high values of ΩGW generated by occasional nearby binaries: the median spectrum always lies within the 1σ CL range. All in all, consistently with earlier Monte Carlo studies of the GW signal (Sesana et al. 2008; Kocsis & Sesana 2011), we find that typical spectra tend to fall below the Ω GW f 2 3 $ \Omega_{\mathrm{GW}} \propto f^{\frac23} $ expectation in higher frequency bins, while a few bins display sharp peaks.

thumbnail Fig. 5.

GW energy spectrum ΩGW and the fractional circular polarization 𝒫circ calculated from 1600 Monte Carlo realizations of the SMBH binary population. In the left panels, the black solid curve shows the mean and the gray bands the 68% and 95% CL regions of the Monte Carlo realizations, and the red curve corresponds to one of the realizations shown in Fig. 3. In the upper left panel the green line shows the mean GW energy density (8). The distributions of ΩGW and 𝒫GW at 1 nHz (blue) and at 30 nHz (orange) calculated from the Monte Carlo realizations are shown in the right panels. In the upper right panel, the dashed lines show P Ω 5 2 $ P \propto \Omega^{-\frac52} $, and in the lower right panel, the dashed curves show Gaussian fits to these distributions of width σ = 0.13 (blue) and σ = 0.23 (orange). In this figure, pBH = 0.17 is assumed.

As there are a finite number of sources, Eq. (8) is subject to statistical fluctuations, which arise mostly from the possibility of having a few strong sources nearby. To obtain an order-of-magnitude estimate, we focus on the closest binaries and ignore the redshift dependence. In this case, Ω GW ( 1 ) D L 2 $ \Omega^{(1)}_{\mathrm{GW}} \propto D_{L}^{-2} $, by Eqs. (8) and (10), and the probability of finding an event at DL is P (1) ( D L ) D L 2 $ P^{(1)}(D_{L}) \propto D_{L}^2 $. Thus Ω GW ( 1 ) D L , min D L , max d D L P ( D L ) Ω GW ( 1 ) D L , max $ \langle \Omega^{(1)}_{\mathrm{GW}} \rangle \propto \int^{D_{L, \rm max}}_{D_{L, \rm min}} {\text{ d}}D_{L} P(D_{L}) \Omega^{(1)}_{\mathrm{GW}} \propto D_{L, \rm max} $, where, in this simplified approach, DL,max is some large luminosity distance at which Ω GW ( 1 ) $ \Omega^{(1)}_{\mathrm{GW}} $ gets suppressed, and DL,min is the distance to the nearest possible source. On the other hand, ( Ω GW ( 1 ) ) 2 1 / D L , min $ \langle (\Omega^{(1)}_{\mathrm{GW}})^{2} \rangle \propto 1/D_{L, \rm min} $. Therefore, we expect the mean of the GW signal to be determined by faraway sources, while the variance is set by a few close-by binaries. Moreover, we can estimate from P(1)(DL) that Ω GW ( 1 ) $ \Omega^{(1)}_{\mathrm{GW}} $ has a relatively flat power-law tail at large values,

P ( 1 ) ( Ω ) Ω 5 2 , when Ω . $$ \begin{aligned} P^{(1)}(\Omega ) \propto \Omega ^{-\frac{5}{2}}, \quad \text{ when} \quad \Omega \rightarrow \infty . \end{aligned} $$(16)

Since the closest distance to massive BH binaries is constrained, this tail will be cut off at Ω GW , max ( 1 ) D L , min 2 $ \Omega^{(1)}_{\mathrm{GW, max}} \propto D_{L, \rm min}^{-2} $. Even with this cutoff, the mean and the variance are still not very useful characteristics of the uncertainties, and we find it more illuminating to estimate the confidence intervals around the median value.

The upper right panel of Fig. 5 displays the distributions of ΩGW at 1 nHz (blue) and at 30 nHz (orange) found in 1600 Monte Carlo realizations of the SMBH binary population. We note that the distributions at both frequencies have tails that approach the analytical result P Ω 5 2 $ P \propto \Omega^{-\frac52} $ (Eq. (16)), as indicated by the dashed lines. These tails exhibit explicitly why the mean (black solid) line in Fig. 5 tends to lie above the 68% CL band, while the median (gray dashed) line lies within it at all frequencies. We note that the overall shapes of the distributions resemble the analytic calculations shown in the left panel of Fig. A.1. In summary, even if a bin has a high number of contributing events (i.e., N ¯ $ \bar N \to \infty $), the distribution of ΩGW does not converge to a Gaussian because it retains its Ω GW 5 / 2 $ \Omega_{\mathrm{GW}}^{-5/2} $ tail. Since the variance is not well behaved, we define the width of the distribution ΔΩGW as the width of the 68% confidence interval, as in Fig. 5. The width-to-mean ratio scales as:

Δ Ω GW / Ω GW N ¯ 1 3 , $$ \begin{aligned} \Delta _{\Omega _{\rm GW}}/\left\langle \Omega _{\rm GW}\right\rangle \propto \bar{N}^{-\frac{1}{3}} , \end{aligned} $$(17)

when N ¯ 1 $ \bar N \gg 1 $, as is shown in Appendix A. This scaling is slower than the typical 1 / N ¯ $ 1/\sqrt{\bar N} $ scaling predicted by the central limit theorem.

3.3. GW polarization

Recent studies have argued that the circular polarization of the signal can be used to estimate whether the SGWB comes from a handful of sources or a relatively large population of binaries (Kato & Soda 2016; Conneely et al. 2019; Hotinli et al. 2019; Belgacem & Kamionkowski 2020; Sato-Polito & Kamionkowski 2022; Valbusa Dall’Armi et al. 2023). We recall that the left and right circular GW polarization amplitudes from a binary with inclination angle θ are

| h L , R ( f ) | 2 = 1 16 ( 1 ± cos θ ) 4 | h ( f ) | 2 , $$ \begin{aligned} |h_{L,R}(f)|^2 = \frac{1}{16}(1 \pm \cos {\theta })^4 \, |h(f)|^2 , \end{aligned} $$(18)

where

| h ( f ) | = 4 M z 5 3 D L ( π f ) 2 3 $$ \begin{aligned} |h(f)| = \frac{4 \mathcal{M} _z^{\frac{5}{3}}}{D_L} (\pi f)^{\frac{2}{3}} \end{aligned} $$(19)

is the maximal GW strain from an inspiralling binary, and the gravitational Stokes parameters are defined by

I ( f ) = | h L ( f ) | 2 + | h R ( f ) | 2 , V ( f ) = | h L ( f ) | 2 | h R ( f ) | 2 . $$ \begin{aligned} \begin{aligned} I(f)&= |h_L(f)|^2 + |h_R(f)|^2 , \\ V(f)&= |h_L(f)|^2 - |h_R(f)|^2 . \end{aligned}\end{aligned} $$(20)

The fractional amount of circular polarization of the SGWB can be characterized by the quantity

P circ ( f ) = i V i ( f ) i I i ( f ) , $$ \begin{aligned} \mathcal{P} _{\rm circ}(f) = \frac{\sum _i V_i(f)}{\sum _i I_i(f)}, \end{aligned} $$(21)

where the sums are over all binaries in a fixed frequency range.

If the GW signal is dominated by a single source, then 𝒫circ depends only on the inclination angle and large circular polarizations are preferred, with 𝒫circ > 0.87(0.2) at the 68% CL (95% CL). On the other hand, if the GW signal is dominated by several (Ndom ≳ 10) sources of comparable strengths, then the 𝒫circ(f) will be approximately Gaussian with a width determined by fluctuations in σV. Since σV/⟨Iθ = 1.17 for a single source, we can estimate that (see Appendix A for details)

σ P circ σ V / I θ 1.17 / N dom . $$ \begin{aligned} \sigma _{\mathcal{P} _{\rm circ}} \approx \sigma _{V}/\langle I\rangle _{\theta } \approx 1.17/\sqrt{N_{\rm dom}} . \end{aligned} $$(22)

The lower-left panel of Fig. 5 illustrates the distributions of the circular polarization 𝒫circ(f) found in the sample of 1600 Monte Carlo realizations of the SMBH binary population with M > 109M and η > 0.01 whose frequency spectra were illustrated in Fig. 5. We see that the mean value of 𝒫circ(f)≈0, as expected, but large statistical fluctuations are possible even at the 1σ level. This phenomenon is visible in the red curve, which shows results from one of the Monte Carlo realizations. The large fluctuations reflect the fact that the observable GW spectrum could be due to a very limited number of sources, particularly at higher frequencies.

The lower-right panel of Fig. 5 displays the distributions of 𝒫circ at 1 nHz (blue) and at 30 nHz (orange) found in 1600 Monte Carlo realizations of the SMBH binary population. We note that the statistical distribution is indeed broader at the higher frequency, as expected. The dashed curves are Gaussian fits to polarization distributions with widths σ = 0.13 (blue) and σ = 0.23 (orange). These fits are very accurate as also demonstrated in the right panel of Fig. A.1.

As the variance of the circular polarization of the signal depends on the number of binaries that dominate the signal, it would be suppressed if the merger rate of the heaviest binaries were suppressed, which could be the case if pBH is not universal. By suppressing pBH, for instance, at M > 109 M and enhancing pBH for 109 M > M > 106 M, we can accommodate the IPTA common-spectrum effect. In this case, the GW energy spectrum ΩGW becomes smoother and approaches the naive Ω GW f 2 3 $ \Omega_{\mathrm{GW}} \propto f^{\frac23} $ behavior while suppressing the circular polarization 𝒫circ. Moreover, in this case, pBH for the lighter binaries would be enhanced, which would increase the number of signals in the sensitivity range of LISA and possibly AEDGE.

3.4. Prospects for future GW observatories

As discussed earlier, BH binaries can generate two qualitatively distinct classes of signal: long, nearly monochromatic signals from the slowly-evolving inspiralling phase and relatively short signals from the merger and ringdown. The detections of these two types of sources need to be considered separately.

The number of detectable nearly monochromatic GW signals that arise from inspiralling BH binaries is

N insp = T + τ min d τ d λ p det [ SNR c SNR ( τ , z ) ] , $$ \begin{aligned} N_{\rm insp} = \int _{\mathcal{T} +\tau _{\rm min}}^\infty \!\!\mathrm{d}\tau \int \mathrm{d}\lambda \,p_{\rm det}\!\left[\frac{\mathrm{SNR}_c}{\mathrm{SNR}(\tau ,z)}\right] , \end{aligned} $$(23)

where we account only for binaries whose coalescence time τ is longer than τmin = 1 day and we estimate the SNR for a detector characterized by the noise power spectrum Sn(f) as

SNR ( τ , z ) = 0 T d t 2 | h ( τ t , z ) | 2 S n ( f ( τ t , z ) ) . $$ \begin{aligned} \mathrm{SNR}(\tau ,z) = \sqrt{\int _0^\mathcal{T} \mathrm{d}t \,\frac{2 |h(\tau -t,z)|^2}{S_n(f(\tau -t,z))}} . \end{aligned} $$(24)

The optimal time-dependent inspiral strain |h(t)| = |h(f(t))| is given by Eq. (19) and the GW frequency as a function of time is given by Eq. (12). Analogously, the expected number of BH binary merger events is

N merg = T d λ p det [ SNR c SNR ( z ) ] , $$ \begin{aligned} N_{\rm merg} = \mathcal{T} \int \mathrm{d}\lambda \,p_{\rm det}\!\left[\frac{\mathrm{SNR}_c}{\mathrm{SNR}(z)}\right] , \end{aligned} $$(25)

where

SNR ( z ) = f ( τ min ) d f 4 | h ( f ) | 2 S n ( f ) . $$ \begin{aligned} \mathrm{SNR}(z) = \sqrt{\int _{f(\tau _{\rm min})}^\infty \mathrm{d}f \,\frac{4 |\tilde{h}(f)|^2}{S_n(f)}}. \end{aligned} $$(26)

The cut at f(τmin) here implies that we only consider the SNR only from the last day of the signal. In both cases, the detection probability pdet accounts for the detector’s antenna patterns and includes the average over the binary inclination, sky location, and polarization (Finn & Chernoff 1993; Gerosa et al. 2019), the noises Sn include the foregrounds from stellar mass BH binaries and white dwarf binaries (see e.g., Lewicki & Vaskonen 2021). Furthermore, we used SNRc = 8 for the detection threshold and 𝒯 = 1 yr for the observation time.

In the monochromatic limit, Eq. (24) is simplified to SNR = | h ( f ) | 2 T / S n ( f ) $ \mathrm{SNR} = |h(f)| \sqrt{2\mathcal{T}/S_n(f)} $. Using this, we can convert the prospected noise Sn(f) of IPTA30 to a lower bound on the strain |h(f)| for which the SNR exceeds the detection threshold. As indicated by the very long coalescence times in the middle right panel of Fig. 3, the signals in the IPTA30 band can be adequately approximated as monochromatic. In the right panel of Fig. 3 the dashed black curve shows the GW strain |h(f)| that for 20 year observation time with IPTA30 gives SNR = 8. We see that IPTA30 can potentially resolve several, 𝒪(10), monochromatic GW signals from SMBH binaries.

At higher frequencies, f > 10−6 Hz, we find that LISA is the only detector capable of measuring the near-monochromatic sources arising in the early stages of the inspiral phase, namely, more than one day from the beginning of the merger phase. It is clear that AEDGE is not able to observe these near-monochromatic sources because the coalescence time of binaries heavier than 103M is less than a day when they enter the AEDGE sensitivity window6. Figure 6 illustrates the prospects to detect the inspiral signals with LISA. The shading in the panels correspond to the number of events whose parameters fall within each rectangular bin. The distribution peaks below ℳ0 = 104 M reflecting the binary population, that increases towards lower masses, and the cut-off at 103M. Similarly, because the GW signal (Eq. (10)) from inspirals does not explicitly depend on η, the η distribution in the right panel of Fig. 6 arises solely from the binary population. The binary merger rate peaks at 4 < z < 6 (see Fig. 1) but LISA cannot spot the majority of the binaries beyond z ∼ 4.

thumbnail Fig. 6.

Expected numbers of near-monochromatic GW signals generated more than a day before the merger that would be detectable in a year of LISA observation, as functions of the chirp mass ℳ0 and redshift z (left panel) and of the chirp mass and symmetric mass ratio η (right panel). The merger rates assume the abundance of IMBHs described by Eq. (6) for BHs heavier than 103M and a universal merger efficiency factor of pBH = 0.17 inferred from IPTA data, as discussed in the text.

We find that, at z < 10, there are a total of 𝒪(107) near-monochromatic sources in the frequency range of 10−5 Hz < f < 0.1 Hz7. As only a small fraction of these can be resolved by LISA, we expect that the rest constitute a significant stochastic GW background. It should be noted that the potential redshift dependence of pBH introduces additional uncertainties in the merger rate at higher redshifts (z ≳ 3), which is not well probed by the PTA measurements. We leave a detailed study of this background and its detectability with LISA for future work.

Figure 7 illustrates AEDGE and LISA prospects for observing GW events from BHs less than a day after the merger. The expected total number of detectable events is slightly larger for LISA than for AEDGE. Due to the different frequency ranges, LISA can spot heavier binaries, ℳ0 > 104M, while AEDGE can see lighter binaries, ℳ0 < 2 × 103M. LISA could (in principle) detect even heavier mergers, namely, of ℳ0 ≲ 108M, but since such mergers are so rare it is unlikely that LISA will see any mergers above 106M. The expected number of detectable mergers peaks at 4 < z < 6 for both of the experiments. For AEDGE, the z distribution of the detectable mergers reflects the binary merger rate (see Fig. 1) while for LISA the z distribution peaks at slightly lower z because the SNR of ℳ0 ∼ 104 binaries for LISA is not as high as it is for AEDGE. The lower panels show that, for all detectors, the BH mass ratio peaks at the highest values. For AION-km and ET we find that the total expected number of events per year is less than three.

thumbnail Fig. 7.

Expected numbers of detectable GW events generated during the last day before the merger in a year of observation by LISA (left panels) and AEDGE (right panels) as functions of the chirp mass ℳ0 and redshift z (upper panels) and of chirp mass and symmetric mass ratio η (lower panels). The merger rates assume the abundance of IMBHs described by Eq. (6) for BHs heavier than 103M and a universal merger efficiency factor of pBH = 0.17 inferred from the IPTA data.

We note that most of the IMBH events detectable by AEDGE during the last day prior to the merger (see the left panels in Fig. 7) will also have been detected by LISA during the previous infall stage, as seen in Fig. 6. This opens up prospects for using LISA data to predict when and in what direction AEDGE will observe IMBH mergers, sharpening tests of general relativity and giving advance warnings for searches for possible multi-messenger signals.

We caution that there are observational uncertainties in the low mass cut-off in Eq. (6), due to the difficulty of measuring very faint active galactic nuclei (AGNs) or inactive BHs in dwarf galaxies. As mentioned above, surveys such as eRASS and AMUSE (Miller et al. 2015) are already constraining this low-mass region and are compatible with the assumed cut-off, although a heavier mass cut-off may be favored (Chadayammuri et al. 2023). This could easily be achieved in models with modified initial fluctuation spectra (Hütsi et al. 2023).

The cut-off mass is tightly related to the SMBH formation mechanism since a BH in the halo centre cannot be lighter than the seed that originated the growth. For a full review, we refer to Volonteri et al. (2021). In this sense, the cut-off we have considered assumes the existence of some light seeds, which augment the possible GW signal. Independently of the growth through halo merging, IMBHs can also be born inside dense stellar media like nuclear and globular star clusters. The growth from stellar masses happens because of repeated encounters inside these dense environments. Estimating the GW emission from such encounters is an active field of research (Fragione & Kocsis 2018; Fragione et al. 2019; Fragione & Loeb 2023) and their signals may contribute significantly to the expected number of IMBH mergers.

4. Conclusions

Here, we describe the results from a model for the PTA nHz common-process signal based on a simulation of massive BH mergers. The magnitude of the signal depends on the merger probability, pBH, which is a product of the probabilities that a pair of halos contain massive BHs and the probability that they will merge. We sidestepped the considerable uncertainties in modeling these probabilities by fitting a mass-independent value of pBH to the PTA signal, finding pBH = 0.17, with a factor of two uncertainty. With this assumption, the dominant contribution to the PTA signal is made by mergers with total masses > 109M and about 10% from masses < 109M. The PTA mergers would have redshifts 𝒪(1) and have mass asymmetries ≳10.

The number of mergers contributing most of the PTA signal at frequencies 𝒪(10) nHz is limited. Consequently, the frequency spectrum becomes quite irregular, the spectral index may deviate from the analytic value of 2/3, individual mergers may be distinguished, and there may be detectable circular polarization. These will be interesting targets for future experiments in the nHz range, including PTAs and SKA.

Assuming the same mass-independent value of pBH as for the PTA signal, there would be observable signals from mergers with total masses ∈(103, 106) M in the LISA experiment and from mergers with total masses ∈(103, 105) M in the AEDGE experiment. Data from these experiments will be able to check the accuracy of the constancy of pBH. In principle, their signals could be even larger than our estimates if pBH is closer to unity for masses below the PTA range, or if there is a significant GW contribution seeded by Population III stars. However, a lower value for pBH and, hence, a lower event rate cannot be excluded.


1

However, alternative scenarios for quasar optical variability may also be considered, such as intrinsic variability in the accretion disc (King et al. 2013).

2

A similar signal could be generated by primordial BHs (Vaskonen & Veermäe 2021; De Luca et al. 2021; Kohri & Terada 2021; Ashoorioon et al. 2022), but might require modifications of models based on simple cosmic string networks (Ellis & Lewicki 2021; Blasi et al. 2021; Buchmuller et al. 2020): see, e.g., Blanco-Pillado et al. (2021).

3

We note in passing that binaries capable of generating the PTA nHz background could not explain the year-like quasi-periodicity seen in blazars (Holgado et al. 2018).

4

In addition to the experiments shown in Fig. 2, we have also considered the prospective sensitivities to GWs using the astrometric data from the Gaia and Nancy Grace Roman space telescopes (Wang et al. 2021). Their nominal sensitivities lie well above the black dashed line in Fig. 2 corresponding to pBH = 1, but a possible improvement of the Roman sensitivity might enable the detection of SMBH binary signals at frequencies ∈(10−7, 10−6) Hz. Moreover, while the prospective power-law integrated sensitivity to GWs through binary resonance by laser ranging of the Moon (Blas & Jenkins 2022) reaches the gray band shown in Fig. 2, the most likely signal at f ∼ 10−6 Hz is well below that (see Sect. 3.2) and consists of a few nearly monochromatic sources.

5

We remark that, in Eq. (11), the observed time spent in a given frequency band is redshifted when compared to the time in the binary rest frame. This redshift has been accounted for in Eq. (12).

6

In this work we have cut the AEDGE sensitivity at 3 mHz. This cut is motivated by the potential Newtonian gravity backgrounds (Hogan et al. 2011; Graham et al. 2017). We note, however, that dedicated studies of these backgrounds for AEDGE have not been performed.

7

We emphasize, however, that there could be lower-mass BHs that are remnants of Population III stars, which are not included in our analysis.

Acknowledgments

This work was supported by European Regional Development Fund through the CoE program grant TK133 and by the Estonian Research Council grants PRG803 and PSG869. The work of J.E. was supported by the United Kingdom STFC Grant ST/T000759/1. The work of V.V. was supported by the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101065736.

References

  1. Ajith, P., Babak, S., Chen, Y., et al. 2008, Phys. Rev. D, 77, 104017 [NASA ADS] [CrossRef] [Google Scholar]
  2. Antoniadis, J., Arzoumanian, Z., Babak, S., et al. 2022, MNRAS, 510, 4873 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arzoumanian, Z., Baker, P. T., Blecha, L., et al. 2023, ApJL, submitted [arXiv:2301.03608] [Google Scholar]
  5. Ashoorioon, A., Rezazadeh, K., & Rostami, A. 2022, Phys. Lett. B, 835, 137542 [NASA ADS] [CrossRef] [Google Scholar]
  6. Badurina, L., Bentine, E., Blas, D., et al. 2020, JCAP, 05, 011 [CrossRef] [Google Scholar]
  7. Badurina, L., Buchmueller, O., Ellis, J., et al. 2021, Phil. Trans. A. Math. Phys. Eng. Sci., 380, 20210060 [Google Scholar]
  8. Barkana, R., & Loeb, A. 2001, Phys. Rept., 349, 125 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022, ApJ, 941, 119 [CrossRef] [Google Scholar]
  10. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  11. Belgacem, E., & Kamionkowski, M. 2020, Phys. Rev. D, 102, 023004 [NASA ADS] [CrossRef] [Google Scholar]
  12. Benson, A. J., Kamionkowski, M., & Hassani, S. H. 2005, MNRAS, 357, 847 [NASA ADS] [CrossRef] [Google Scholar]
  13. Berti, E., Cardoso, V., & Will, C. M. 2006, Phys. Rev. D, 73, 064030 [NASA ADS] [CrossRef] [Google Scholar]
  14. Blanco-Pillado, J. J., Olum, K. D., & Wachter, J. M. 2021, Phys. Rev. D, 103, 103512 [NASA ADS] [CrossRef] [Google Scholar]
  15. Blas, D., & Jenkins, A. C. 2022, Phys. Rev. Lett., 128, 101103 [NASA ADS] [CrossRef] [Google Scholar]
  16. Blasi, S., Brdar, V., & Schmitz, K. 2021, Phys. Rev. Lett., 126, 041305 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440 [NASA ADS] [CrossRef] [Google Scholar]
  18. Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53 [Google Scholar]
  19. Buchmuller, W., Domcke, V., & Schmitz, K. 2020, Phys. Lett. B, 811, 135914 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chadayammuri, U., Bogdan, A., Ricarte, A., & Natarajan, P. 2023, ApJ, 946, 51 [CrossRef] [Google Scholar]
  21. Chen, S., Middleton, H., Sesana, A., Del Pozzo, W., & Vecchio, A. 2017, MNRAS, 468, 404 [NASA ADS] [CrossRef] [Google Scholar]
  22. Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970 [NASA ADS] [CrossRef] [Google Scholar]
  23. Conneely, C., Jaffe, A. H., & Mingarelli, C. M. F. 2019, MNRAS, 487, 562 [NASA ADS] [CrossRef] [Google Scholar]
  24. De Luca, V., Franciolini, G., & Riotto, A. 2021, Phys. Rev. Lett., 126, 041303 [NASA ADS] [CrossRef] [Google Scholar]
  25. De Paolis, F., Ingrosso, G., & Nucita, A. A. 2002, A&A, 388, 470 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. De Paolis, F., Ingrosso, G., & Nucita, A. A. 2004, A&A, 426, 379 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Dodelson, S. 2003, Modern Cosmology (Amsterdam: Academic Press) [Google Scholar]
  28. Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [Google Scholar]
  29. El-Neaj, Y. A., Alpigiani, C., Amairi-Pyka, S., et al. 2020, EPJ Quant. Technol., 7, 6 [CrossRef] [Google Scholar]
  30. Ellis, J., & Lewicki, M. 2021, Phys. Rev. Lett., 126, 041304 [NASA ADS] [CrossRef] [Google Scholar]
  31. Enoki, M., & Nagashima, M. 2007, Prog. Theor. Phys., 117, 241 [CrossRef] [Google Scholar]
  32. Enoki, M., Inoue, K. T., Nagashima, M., & Sugiyama, N. 2004, ApJ, 615, 19 [NASA ADS] [CrossRef] [Google Scholar]
  33. Erickcek, A. L., Kamionkowski, M., & Benson, A. J. 2006, MNRAS, 371, 1992 [NASA ADS] [CrossRef] [Google Scholar]
  34. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
  35. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
  36. Finn, L. S., & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fragione, G., & Kocsis, B. 2018, Phys. Rev. Lett., 121, 161103 [Google Scholar]
  38. Fragione, G., & Loeb, A. 2023, ApJ, 944, 81 [CrossRef] [Google Scholar]
  39. Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47 [NASA ADS] [CrossRef] [Google Scholar]
  40. Gerosa, D., Ma, S., Wong, K. W. K., et al. 2019, Phys. Rev. D, 99, 103004 [NASA ADS] [CrossRef] [Google Scholar]
  41. Goncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021, ApJ, 917, L19 [NASA ADS] [CrossRef] [Google Scholar]
  42. Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74 [Google Scholar]
  43. Graham, P. W., Hogan, J. M., Kasevich, M. A., Rajendran, S., & Romani, R. W. 2017, ArXiv e-prints [arXiv:1711.02225] [Google Scholar]
  44. Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
  45. Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hogan, J. M., Johnson, D. M. S., Dickerson, S., et al. 2011, Gen. Rel. Grav., 43, 1953 [NASA ADS] [CrossRef] [Google Scholar]
  47. Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, MNRAS, 481, L74 [NASA ADS] [CrossRef] [Google Scholar]
  48. Hotinli, S. C., Kamionkowski, M., & Jaffe, A. H. 2019, Open J. Astrophys., 2, 8 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hütsi, G., Raidal, M., Urrutia, J., Vaskonen, V., & Veermäe, H. 2023, Phys. Rev. D, 107, 043502 [CrossRef] [Google Scholar]
  50. Iguchi, S., Okuda, T., & Sudou, H. 2010, ApJ, 724, L166 [NASA ADS] [CrossRef] [Google Scholar]
  51. Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2021, MNRAS, 509, 3488 [CrossRef] [Google Scholar]
  52. Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, PoS, AASKA14, 037 [Google Scholar]
  53. Kaiser, A. R., & McWilliams, S. T. 2021, CQG, 38, 055009 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kato, R., & Soda, J. 2016, Phys. Rev. D, 93, 062003 [NASA ADS] [CrossRef] [Google Scholar]
  55. Kelley, L. Z., Blecha, L., & Hernquist, L. 2017a, MNRAS, 464, 3131 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017b, MNRAS, 471, 4508 [NASA ADS] [CrossRef] [Google Scholar]
  57. King, O. G., Hovatta, T., Max-Moerbeck, W., et al. 2013, MNRAS, 436, 114 [Google Scholar]
  58. Kocsis, B., & Sesana, A. 2011, MNRAS, 411, 1467 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kohri, K., & Terada, T. 2021, Phys. Lett. B, 813, 136040 [NASA ADS] [CrossRef] [Google Scholar]
  60. Kovačević, A. B., Songsheng, Y.-Y., Wang, J.-M., & Popovic, L. C. 2022, A&A, 663, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Lacey, C. G., & Cole, S. 1993, MNRAS, 262, 627 [NASA ADS] [CrossRef] [Google Scholar]
  62. Lewicki, M., & Vaskonen, V. 2021, Eur. Phys. J. C, 83, 168 [Google Scholar]
  63. LIGO Scientific & Virgo Collaborations (Abbott, B. P., et al.) 2016, Phys. Rev. Lett., 116, 061102 [CrossRef] [PubMed] [Google Scholar]
  64. LIGO Scientific, Virgo & KAGRA Collaborations (Abbott, R., et al.) 2021, ArXiv e-prints [arXiv:2111.03634] [Google Scholar]
  65. Lippai, Z., Frei, Z., & Haiman, Z. 2009, ApJ, 701, 360 [CrossRef] [Google Scholar]
  66. LISA Collaboration (Amaro-Seoane, P., et al.) 2017, ArXiv e-prints [arXiv:1702.00786] [Google Scholar]
  67. Merritt, D. 2013, CQG, 30, 244005 [NASA ADS] [CrossRef] [Google Scholar]
  68. Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98 [NASA ADS] [CrossRef] [Google Scholar]
  69. O’Neill, S., Kiehlmann, S., Readhead, A. C. S., et al. 2022, ApJ, 926, L35 [CrossRef] [Google Scholar]
  70. Peters, P. C. 1964, Phys. Rev., 136, B1224 [NASA ADS] [CrossRef] [Google Scholar]
  71. Peters, P. C., & Mathews, J. 1963, Phys. Rev., 131, 435 [Google Scholar]
  72. Phinney, E. S. 2001, ArXiv e-prints [arXiv:astro-ph/0108028] [Google Scholar]
  73. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
  75. Rieger, F. M., & Mannheim, K. 2000, A&A, 359, 948 [NASA ADS] [Google Scholar]
  76. Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, CQG, 29, 124013 [NASA ADS] [CrossRef] [Google Scholar]
  77. Sato-Polito, G., & Kamionkowski, M. 2022, Phys. Rev. D, 106, 023004 [NASA ADS] [CrossRef] [Google Scholar]
  78. Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623 [NASA ADS] [CrossRef] [Google Scholar]
  79. Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192 [NASA ADS] [CrossRef] [Google Scholar]
  80. Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255 [NASA ADS] [CrossRef] [Google Scholar]
  81. Taylor, S. R., Huerta, E. A., Gair, J. R., & McWilliams, S. T. 2016, ApJ, 817, 70 [NASA ADS] [CrossRef] [Google Scholar]
  82. Valbusa Dall’Armi, L., Nishizawa, A., Ricciardone, A., & Matarrese, S. 2023, Phys. Rev. Lett., 131, 041401 [CrossRef] [Google Scholar]
  83. Valtonen, M. J., et al. 2008, Nature, 452, 851 [NASA ADS] [CrossRef] [Google Scholar]
  84. Vaskonen, V., & Veermäe, H. 2021, Phys. Rev. Lett., 126, 051303 [NASA ADS] [CrossRef] [Google Scholar]
  85. Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [Google Scholar]
  86. Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nat. Rev. Phys., 3, 732 [NASA ADS] [CrossRef] [Google Scholar]
  87. Wang, Y., Pardo, K., Chang, T.-C., & Doré, O. 2021, Phys. Rev. D, 103, 084007 [NASA ADS] [CrossRef] [Google Scholar]
  88. Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 590, 691 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Statistics of the GW energy spectrum and circular polarization

A.1. Distribution of ΩGW(f)

The GW energy spectrum arising from a population of BH binaries can be expressed as the sum of the contributions from individual binaries. We first consider the contribution to ΩGW coming from a single inspiralling binary:

Ω GW ( 1 ) ( θ ) = | d τ d ln f | 1 2 π 5 f 3 | h ( f ) | 2 ρ c , $$ \begin{aligned} \Omega _{\rm GW}^{(1)}(\boldsymbol{\theta }) = \left| \frac{\mathrm{d}\tau }{\mathrm{d}\ln f} \right|^{-1} \frac{2\pi }{5} \frac{f^3 |\tilde{h}(f)|^2}{\rho _c} \, , \end{aligned} $$(A.1)

where θ ≡ {ℳz,z,f} are parameters describing the binary. The GW energy spectrum arising from the inspiralling binary population (or any of its subpopulations) arises from the sum of nearly monochromatic components

Ω GW ( f ) i Ω GW ( 1 ) ( θ i ) δ ( ln f / f i ) . $$ \begin{aligned} \Omega _{\rm GW}(f) \approx \sum _{i} \Omega _{\mathrm{GW}}^{(1)}(\boldsymbol{\theta }_i) \delta (\ln f/f_i) \, . \end{aligned} $$(A.2)

This approximation is valid as long as the change in frequency over the observation period is smaller than the spectral resolution. The parameters θi as well as Ω GW ( 1 ) $ \Omega_{\mathrm{GW}}^{(1)} $ are independent and identically distributed. Thus, the statistical properties of ΩGW(f) can be inferred from the distribution of Ω GW ( 1 ) ( f ) $ \Omega_{\mathrm{GW}}^{(1)}(f) $, which is given by

P ( 1 ) ( Ω | f ) = 1 n ( f ) d λ | d τ d ln f | δ ( Ω Ω GW ( 1 ) ) = F ( Ω f 10 3 ) f Ω 3 2 n ( f ) , $$ \begin{aligned} \begin{aligned} P^{(1)}(\Omega |f)&= \frac{1}{n(f)} \int \mathrm{d}\lambda \left| \frac{\mathrm{d}\tau }{\mathrm{d}\ln f} \right| \delta \left( \Omega - \Omega ^{(1)}_{\rm GW} \right) \\&= \frac{\mathcal{F} (\Omega f^{-\frac{10}{3}})}{f \Omega ^{\frac{3}{2}} n(f) } \, , \end{aligned}\end{aligned} $$(A.3)

where we define the following function:

F ( x ) d z d η [ λ ( M z , z ) M z η z × M z 16 10 ρ c D L ] M z = 1 π ( 5 π 8 D L 3 ρ c x ) 3 10 $$ \begin{aligned} \begin{aligned} \mathcal{F} (x) \equiv \int \mathrm{d}z&\mathrm{d}\eta \bigg [ \frac{\partial \lambda (\mathcal{M} _z,z)}{\partial \mathcal{M} _z \partial \eta \partial z} \\&\ \times \frac{\mathcal{M} _z}{16 \sqrt{10 \rho _c} D_L} \, \bigg ]_{\mathcal{M} _z = \frac{1}{\pi } \left(\frac{5\pi }{8}D_L^3 \rho _c x \right)^{\frac{3}{10}}} \end{aligned}\end{aligned} $$(A.4)

and the spectral source density:

n ( f ) d M z d η d z d N d M z d η d z d ln f . $$ \begin{aligned} n(f) \equiv \int \mathrm{d}\mathcal{M} _z \mathrm{d}\eta \mathrm{d}z \frac{\mathrm{d}N}{\mathrm{d}\mathcal{M} _z \mathrm{d}\eta \mathrm{d}z \mathrm{d}\ln f} \,. \end{aligned} $$(A.5)

We see that, for circular inspiralling binaries, the four-parameter distribution (11) can be reduced to two functions of a single parameter. Further reductions are unlikely, as these functions depend on the model of the binary merger rate, which we assume to stem from the merger rate of galaxies.

The average energy spectrum in a frequency bin (fj, fj + 1) can be obtained from

Ω GW ( f j ) 1 ln ( f j + 1 / f j ) f j f j + 1 Ω GW ( θ ; f ) d ln f = 1 ln ( f j + 1 / f j ) i = 1 N ( f j ) Ω GW , i ( 1 ) , $$ \begin{aligned} \begin{aligned} \Omega _{\mathrm{GW}}(f_j)&\equiv \frac{1}{\ln (f_{j+1}/f_j)} \int ^{f_{j+1}}_{f_j} \Omega _{\mathrm{GW}}(\boldsymbol{\theta };f) \mathrm{d}\ln f\\&= \frac{1}{\ln (f_{j+1}/f_j)} \sum ^{N(f_j)}_{i=1} \Omega ^{(1)}_{\mathrm{GW},i} \, , \end{aligned}\end{aligned} $$(A.6)

where the number of binaries N(fj) is drawn from a Poisson distribution with the expected value,

N ¯ ( f j ) = f j f j + 1 d ln f n ( f ) , $$ \begin{aligned} \bar{N}(f_j) = \int _{f_{j}}^{f_{j+1}} \!\mathrm{d}\ln f \, n(f) \, , \end{aligned} $$(A.7)

and the Ω GW , i ( 1 ) $ \Omega^{(1)}_{\mathrm{GW},i} $ are drawn from a random distribution P ( 1 ) ( Ω , f ) = P ( 1 ) ( Ω | f ) n ( f ) / N ¯ ( f j ) $ P^{(1)}(\Omega,f) = P^{(1)}(\Omega|f) n(f)/\bar N(f_j) $ with f ∈ (fj, fj + 1). In the following, we suppress fj.

The moment generating function of ΩGW is given by:

M Ω GW ( s ) exp ( s Ω GW ) = N 0 p N i = 1 N exp ( s Ω GW ( 1 ) ) = N 0 N ¯ N N ! e N ¯ ( M Ω GW ( 1 ) ( s ) ) N ¯ = exp [ N ¯ ( M Ω GW ( 1 ) ( k ) 1 ) ] . $$ \begin{aligned} \begin{aligned} M_{\Omega _{\rm GW}}(s)&\equiv \langle \exp \left(s \Omega _{\rm GW} \right) \rangle \\&= \sum _{N \ge 0} p_N \prod ^{N}_{i=1} \left\langle \exp \left(s\Omega ^{(1)}_{\rm GW} \right) \right\rangle \\&= \sum _{N \ge 0} \frac{{\bar{N}}^{N}}{N!} e^{-\bar{N}}\left( M_{\Omega ^{(1)}_{\rm GW}}(s) \right)^{\bar{N}} \\&= \exp \left[ \bar{N} \left(M_{\Omega ^{(1)}_{\rm GW}}(k) - 1 \right) \right] \, . \end{aligned}\end{aligned} $$(A.8)

Thus, the cumulant generating function of ΩGW(fj), that is,

K Ω GW ( s ) ln M Ω GW ( s ) = f j f j + 1 d N [ exp ( s Ω GW ( 1 ) ) 1 ] . $$ \begin{aligned} \begin{aligned} K_{\Omega _{\rm GW}}(s)&\equiv \ln M_{\Omega _{\rm GW}}(s) \, \\&= \int _{f_{j}}^{f_{j+1}} \mathrm{d}N \left[ \exp (s \,\Omega ^{(1)}_{\rm GW}) - 1 \right] \, . \end{aligned}\end{aligned} $$(A.9)

is proportional to the generating function of M Ω GW ( 1 ) ( s ) $ M_{\Omega^{(1)}_{\mathrm{GW}}}(s) $, implying that the nth cumulant of ΩGW is

κ n [ Ω GW ] = N ¯ ( Ω GW ( 1 ) ) n = f j f j + 1 d N ( Ω GW ( 1 ) ) n . $$ \begin{aligned} \begin{aligned} \kappa _{n} [\Omega _{\rm GW}] = \bar{N} \langle (\Omega ^{(1)}_{\rm GW})^{n} \rangle \, = \int _{f_{j}}^{f_{j+1}} \mathrm{d}N (\Omega ^{(1)}_{\rm GW})^{n} \, . \end{aligned}\end{aligned} $$(A.10)

Therefore, the mean κ1GW)≡⟨ΩGW⟩ matches the expectation value given in Eq. (8), while the variance,

κ 2 ( Ω GW ) ( δ Ω GW ) 2 = N ¯ ( Ω GW ( 1 ) ) 2 $$ \begin{aligned} \kappa _{2} (\Omega _{\rm GW}) \equiv \langle (\delta \Omega _{\rm GW})^{2} \rangle = \bar{N} \langle (\Omega ^{(1)}_{\rm GW})^{2} \rangle \end{aligned} $$(A.11)

is divergent due to the long P ( 1 ) ( Ω ) Ω 5 2 $ P^{(1)}(\Omega) \propto \Omega^{-\frac52} $ tail (16), unless a minimal distance to the closest BH binary is imposed. In particular, as we will demonstrate briefly, this long tail is preserved in the distribution of ΩGW, i.e., the higher cumulants are not diminished when summing several instances of P ( 1 ) ( Ω ) Ω 5 2 $ P^{(1)}(\Omega) \propto \Omega^{-\frac52} $.

The moment-generating function MΩGW(−s) is a Laplace transform of the probability distribution. Thus as an alternative to the Monte Carlo approach adopted in the main text, the probability distribution of ΩGW(fj) can be obtained by an inverse Laplace transform of the moment generating function given in Eq. (A.8). It can be expressed as

P ( Ω GW ) = 1 2 π N ¯ d s e i s Ω GW / N ¯ × × exp ( N ¯ exp ( i s Ω GW ( 1 ) / N ¯ ) 1 ) , $$ \begin{aligned} \begin{aligned} P(\Omega _{\rm GW})&= \frac{1}{2\pi \bar{N}} \int ^{\infty }_{-\infty } \mathrm{d}s \, e^{i s \Omega _{\rm GW}/\bar{N} } \times \\&\times \exp \left( \bar{N} \left\langle \exp \left(-i s \Omega ^{(1)}_{\rm GW}/\bar{N} \right) - 1 \right\rangle \right) \, , \end{aligned}\end{aligned} $$(A.12)

where the average is taken over the single event distribution P(1).

A.1.1. The large N limit. –

The central limit theorem (and Eq. (A.11)) dictates that when N → ∞, then the distribution (A.6) should approach a Gaussian with its relative width scaling as σ Ω GW / Ω GW 1 / N $ \sigma_{\Omega_{\mathrm{GW}}}/\langle \Omega_{\mathrm{GW}} \rangle \propto 1/\sqrt{N} $. However, this is not the case when the second cumulant diverges. Considering a generic distribution with a tail (droping the "GW" subindex for the sake of brevity):

P ( 1 ) ( Ω ) Ω C Ω 5 2 , $$ \begin{aligned} P^{(1)} (\Omega ) \mathop {\sim }\limits ^{\Omega \rightarrow \infty } C \Omega ^{-\frac{5}{2}}, \end{aligned} $$(A.13)

where C is a constant. The second moment ⟨Ω2⟩ diverges. We assume that the small Ω behavior is such that the first moment Ω ¯ ( 1 ) Ω ( 1 ) $ \bar \Omega^{(1)} \equiv \langle \Omega^{(1)}\rangle $ is finite. The N ¯ $ \bar N \to \infty $ asymptotic of the average in Eq. (A.12) can then be computed by separating it into the mean and a term for which the expectation value can be determined by approximating the distribution by its tail (A.13):

N ¯ exp ( i s Ω ( 1 ) N ¯ ) 1 = i s Ω ¯ ( 1 ) + N ¯ exp ( i s Ω ( 1 ) N ¯ ) 1 + i s Ω ( 1 ) N ¯ N i s Ω ¯ ( 1 ) + e i π 3 4 π 4 C 3 N ¯ s 3 2 . $$ \begin{aligned}&\bar{N} \left\langle \exp \left(-i s \frac{\Omega ^{(1)}}{\bar{N}} \right) - 1 \right\rangle \nonumber \\ =&-i s \bar{\Omega }^{(1)} + \bar{N} \left\langle \exp \left(-i s \frac{\Omega ^{(1)}}{\bar{N}} \right) - 1 + i s \frac{\Omega ^{(1)}}{\bar{N}} \right\rangle \nonumber \\ \mathop {\sim }\limits ^{N \rightarrow \infty }&-i s \bar{\Omega }^{(1)} + e^{i \pi \frac{3}{4}} \frac{\sqrt{\pi }4C}{3\sqrt{\bar{N}}} s^{\frac{3}{2}} \, . \end{aligned} $$(A.14)

In summary, the N → ∞ distribution asymptotes to

P ( Ω ) ( C N ¯ ) 2 3 P ( ( C N ¯ ) 2 3 ( Ω Ω ¯ ) ) , $$ \begin{aligned} \begin{aligned} P(\Omega )&\approx (C\bar{N})^{-\frac{2}{3}} \mathcal{P} \left( (C\bar{N})^{-\frac{2}{3}} ( \Omega - \bar{\Omega })\right) \, , \end{aligned}\end{aligned} $$(A.15)

where Ω ¯ N ¯ Ω ¯ ( 1 ) $ \bar{\Omega} \equiv \bar{N}\bar{\Omega}^{(1)} $ denotes the expectation value of Ω and

P ( x ) Re e i π / 6 π 0 d s exp ( i e i π / 6 x s 4 π 3 s 3 2 ) $$ \begin{aligned} \begin{aligned} \mathcal{P} \left( x \right) \equiv \mathrm{Re} \frac{e^{i\pi /6}}{\pi } \int ^{\infty }_{0} \mathrm{d}s \, \exp \left(i e^{i\pi /6} x s - \frac{4\sqrt{\pi }}{3} s^{\frac{3}{2}} \right) \, \end{aligned}\end{aligned} $$(A.16)

is a universal function that is independent of the single event distribution – all information about the latter enters through Ω ¯ ( 1 ) $ \bar{\Omega}^{(1)} $ and C. This limiting case reveals a few crucial features:

  1. Since 𝒫(x)∼x−5/2 when x → ∞, the power-law behavior of the tail is preserved in the large N ¯ $ \bar N $ limit:

    P ( Ω ) Ω Ω ¯ C N ¯ ( Ω Ω ¯ ) 5 2 . $$ \begin{aligned} P(\Omega ) \mathop {\approx }\limits ^{\Omega \gg \bar{\Omega }} \frac{C\bar{N}}{(\Omega - \bar{\Omega })^{\frac{5}{2}}}\, . \end{aligned} $$(A.17)

  2. When the number of contributing events N is increased by a factor of A to AN, the distribution evolves as

    P ( A N ) ( Ω ) = A 2 3 P ( N ) ( A 2 3 ( Ω Ω ¯ ( A N ) ) + Ω ¯ ( N ) ) . $$ \begin{aligned} P^{(AN)}(\Omega )\!=\! A^{-\frac{2}{3}} P^{(N)}\left(\!A^{-\frac{2}{3}}(\Omega \!-\! \bar{\Omega }^{(AN)})\!+\!\bar{\Omega }^{(N)}\!\right). \end{aligned} $$(A.18)

  3. Although the variance diverges, we can define the width of the distribution ΔΩ ≡ Ω+ − Ω from the confidence interval (Ω, Ω+) centred around, e.g., the median. Eq. (A.18) then implies the scaling

    Δ Ω N ¯ 2 / 3 . $$ \begin{aligned} \Delta _{\Omega } \propto \bar{N}^{2/3} \,. \end{aligned} $$(A.19)

    So, the relative width Δ Ω / Ω ¯ $ \Delta_{\Omega}/\bar{\Omega} $ will approach zero, but at a slower pace than when the central limit theorem applies.

  4. Using the method of steepest decent, we find that

    P ( x ) Ω Ω ¯ Ω ¯ Ω 2 π C N ¯ e ( Ω ¯ Ω ) 3 12 π ( C N ¯ ) 2 , $$ \begin{aligned} \mathcal{P} \left( x \right) \mathop {\approx }\limits ^{\Omega \ll \bar{\Omega }} \frac{ \sqrt{\bar{\Omega }-\Omega }}{2 \pi C \bar{N}}e^{-\frac{(\bar{\Omega }-\Omega )^3}{12 \pi (C \bar{N})^2}}\, , \end{aligned} $$(A.20)

    implying values smaller than the mean, that is Ω < Ω ¯ $ \Omega < \bar \Omega $, are much less likely than for Gaussian distributions, as they are suppressed by an exponent of a cube.

As an illustration, we consider a population of sources for which individual signals follow the simple long-tailed distribution

P ( 1 ) ( Ω ) Ω 5 2 θ ( Ω Ω min ) , $$ \begin{aligned} P^{(1)}(\Omega ) \propto \Omega ^{-\frac{5}{2}} \theta (\Omega - \Omega _{\min })\, , \end{aligned} $$(A.21)

which mimics the statistics of ΩGW from a realistic BH population found in section 3.2. In Fig. A.1, we show the resulting distributions of Ω / Ω ¯ $ \Omega/\bar\Omega $ if the signal consists of 10, 100, and 105 sources. The mean is given by Ω ¯ 3 N Ω min $ \bar{\Omega} \equiv 3N\Omega_{\mathrm{min}} $. The distribution is computed by generating random realizations of the source populations (shown by points) and by using (A.12) (shown by solid lines). These two approaches are in excellent agreement. As expected, one can observe that the total signal will inherit the long Ω−5/2 tail from the distribution of individual sources (shown by the dashed line). At large Ω, the tail approaches Eq. (A.17) in all cases, while the asymptotic (A.15) works well for N = 105, but predicts a slightly too wide peak for N ≤ 105. As predicted, the peak gets narrower and moves closer to the expectation value, when N is increased. We note, however, that although the shape of the simplified distribution in Fig. A.1 resembles the distribution in Fig. 5 obtained from a model of the heavy BH binary population, the shape around the peak in the latter depends also on the z dependence and the binary mass distribution.

thumbnail Fig. A.1.

Distributions of the total signal strengths from 1, 10, 100, and 105 sources whose independent signal strengths have identical Ω 5 2 $ \Omega^{-\frac52} $ tails (left). Distribution of 𝒫circ for 1, 10 and 100 identical sources (right). The dashed line shows the distribution from a single source, the solid lines show the analytic estimates, and the distributions obtained from an explicit Monte Carlo analysis are shown by the points.

A.2. Distribution of 𝒫circ

Considering now the distribution of the circular polarization

P circ i = 1 N V i i = 1 N I i , $$ \begin{aligned} \mathcal{P} _{\rm circ} \equiv \frac{\sum ^{N}_{i=1} V_i}{\sum ^{N}_{i=1} I_i} \, , \end{aligned} $$(A.22)

which depends on the distribution and number N of sources. We denote the contributions from individual sources by Vi, Ii given by Eqs. (18-20). As above, we are working in a fixed frequency bin and, thus, suppressing the frequency dependence in the arguments of the functions.

The polarization will vanish on average, that is, at ⟨𝒫circ⟩ = 0, but polarization measurements can provide information from the size of the deviation from that expectation. For a single dominant source, the contribution from the amplitude drops out and we find

P circ = 4 c θ ( 1 + c θ 2 ) 1 + 6 c θ 2 + c θ 4 , $$ \begin{aligned} \mathcal{P} _{\rm circ} = \frac{4 c_{\theta }(1+c_{\theta }^2)}{1 + 6c_{\theta }^2 + c_{\theta }^4} \, , \end{aligned} $$(A.23)

where cθ ≡ cos(θ). The corresponding probability distribution is shown as the dashed curve in Fig. A.1 and can be seen to be dominated by the region |𝒫circ|≈1.

If the signal is dominated by several sources, the analytic treatment of 𝒫circ is complicated due to correlations between V and I. The mean and variance of I and V are given by

I = 4 5 | h | 2 , V = 0 , δ I 2 = 284 315 ( δ | h | 2 ) 2 , δ V 2 = 92 105 ( δ | h | 2 ) 2 , $$ \begin{aligned} \begin{aligned} \langle I \rangle&= \frac{4}{5}\langle |h|^2 \rangle \,, \qquad \qquad \ \langle V \rangle = 0 \,, \\ \langle \delta I^2 \rangle&= \frac{284}{315}\langle (\delta |h|^2)^2 \rangle \,, \ \ \langle \delta V^2 \rangle = \frac{92}{105} \langle (\delta |h|^2)^2 \rangle \,, \end{aligned}\end{aligned} $$(A.24)

where we used the fact that the inclination is independent of other parameters of the binary and the fluctuations in |h|2 are computed as in (A.11). The latter is expected to be large. To estimate the fluctuations in 𝒫circ in the limit when the number of sources is large, we first assume that I can be replaced by its average over inclinations ⟨Iθ ≡ (4/5)∑i|hi|2, but not over signal strengths, so that

P circ V / I θ i = 1 N p i f V ( c θ , i ) , $$ \begin{aligned} \mathcal{P} _{\rm circ} \approx V/\langle I \rangle _{\theta } \equiv \sum ^{N}_{i=1} p_i f_V(c_{\theta ,i}) \,, \end{aligned} $$(A.25)

where pi ≡ (4/5)|hi|2/⟨Iθ ∈ [0, 1] is the fractional contribution to the signal from source i and f V ( c θ )=(5/4) c θ (1+ c θ 2 ) $ f_V(c_{\theta}) = (5/4) c_{\theta}(1+c_{\theta}^2) $ contains the inclination dependence of V. In this way, it is possible to contain the effect of the large fluctuations in ⟨Iθ in the random fractions pi, which vary only in the range [0, 1]. We find:

δ P circ 2 = 115 84 N | p i | 2 . $$ \begin{aligned} \langle \delta \mathcal{P} _{\rm circ}^2 \rangle = \frac{115}{84} N \langle |p_i|^2 \rangle \,. \end{aligned} $$(A.26)

We note that N⟨|pi|2⟩ takes values in the range [1/N, 1] and could be interpreted as a measure of the effective number of dominant sources 1/Ndom. For instance, in the limiting case where the signal is sourced by N binaries that contribute equally, that is, pi = 1/N, we would obtain

σ P circ = 1.17 / N , $$ \begin{aligned} \sigma _{\mathcal{P} _{\rm circ}} = 1.17/\sqrt{N} \, , \end{aligned} $$(A.27)

as also stated in Eq. (22). In the right panel of Fig. A.1, approximate Gaussian distributions with a width given by (A.27) and shown as solid lines are compared to distributions obtained from explicitly generated random populations of 𝒫circ arising from exactly 10 and 100 sources of equal strength (shown by points). We see an excellent agreement between the naive Gaussian approximation and the Monte Carlo estimate already for N = 10.

Finally, we remark that the approximation (A.25) above can be improved by expanding 𝒫circ in δI ≡ I − ⟨Iθ, that is, 𝒫circ = V/⟨Iθ(1 − δI/⟨Iθ + …). This allows us to account for correlations between V and δI, for instance, when computing the variance (A.26),

δ P circ 2 = 115 84 N | p i | 2 2855 1386 N | p i | 3 + $$ \begin{aligned} \langle \delta \mathcal{P} _{\rm circ}^2 \rangle = \frac{115}{84} N \langle |p_i|^2 \rangle - \frac{2855}{1386} N \langle |p_i|^3 \rangle + \ldots \end{aligned} $$(A.28)

The corrections arising from higher powers of δI/⟨Iθ are suppressed by increasing powers of 1/N.

All Figures

thumbnail Fig. 1.

BH merger rates calculated as functions of the redshift z in three ranges of the total merging BH mass M and normalized relative to the merger probability pBH discussed in the text.

In the text
thumbnail Fig. 2.

Mean GW energy density spectrum from massive BH mergers compared with the sensitivities of different experiments. The black dashed curve shows the case where all mergers of galaxies produce a GW signal for BHs heavier than 103M (pBH = 1). The colored bands show the spectra from SMBHs heavier than 109M (dark blue), from SMBHs in the range (106M, 109M); light blue), and from IMBHs with masses in the range 103 − 106M (green), assuming a universal efficiency factor p BH = 0 . 17 0.08 + 0.18 $ p_{\mathrm{BH}} = 0.17^{+0.18}_{-0.08} $. The shaded regions show the prospective sensitivities of IPTA30 (Chen et al. 2017; Kaiser & McWilliams 2021), LISA (LISA Collaboration 2017), AEDGE (El-Neaj et al. 2020; Badurina et al. 2021), AION-km (Badurina et al. 2020, 2021) and ET (Sathyaprakash et al. 2012).

In the text
thumbnail Fig. 3.

Distributions of the binary parameters for a sampling of the BH binary population with 0 < z < 3, 1 nHz < f < 30nHz and M 0 > 10 9 ( f / nHz ) 4 3 $ \mathcal{M}_0 > 10^9 (f/\mathrm{nHz})^{-\frac43} $. The symmetric mass ratio is integrated with the lower bound η > 0.001. In the right panel, the black curve shows the NANOGrav 95% CL constraint on continuous GW sources (Arzoumanian et al. 2023) and the projected sensitivity of IPTA30 at SNR = 8.

In the text
thumbnail Fig. 4.

Expected number of SMBH binaries heavier than Mmin and with symmetric mass ratio η > 0.01 emitting GWs in the indicated frequency range for pBH = 0.17. The vertical dashed lines indicate the fractions of the total GW signal in each frequency band generated by binaries with larger masses.

In the text
thumbnail Fig. 5.

GW energy spectrum ΩGW and the fractional circular polarization 𝒫circ calculated from 1600 Monte Carlo realizations of the SMBH binary population. In the left panels, the black solid curve shows the mean and the gray bands the 68% and 95% CL regions of the Monte Carlo realizations, and the red curve corresponds to one of the realizations shown in Fig. 3. In the upper left panel the green line shows the mean GW energy density (8). The distributions of ΩGW and 𝒫GW at 1 nHz (blue) and at 30 nHz (orange) calculated from the Monte Carlo realizations are shown in the right panels. In the upper right panel, the dashed lines show P Ω 5 2 $ P \propto \Omega^{-\frac52} $, and in the lower right panel, the dashed curves show Gaussian fits to these distributions of width σ = 0.13 (blue) and σ = 0.23 (orange). In this figure, pBH = 0.17 is assumed.

In the text
thumbnail Fig. 6.

Expected numbers of near-monochromatic GW signals generated more than a day before the merger that would be detectable in a year of LISA observation, as functions of the chirp mass ℳ0 and redshift z (left panel) and of the chirp mass and symmetric mass ratio η (right panel). The merger rates assume the abundance of IMBHs described by Eq. (6) for BHs heavier than 103M and a universal merger efficiency factor of pBH = 0.17 inferred from IPTA data, as discussed in the text.

In the text
thumbnail Fig. 7.

Expected numbers of detectable GW events generated during the last day before the merger in a year of observation by LISA (left panels) and AEDGE (right panels) as functions of the chirp mass ℳ0 and redshift z (upper panels) and of chirp mass and symmetric mass ratio η (lower panels). The merger rates assume the abundance of IMBHs described by Eq. (6) for BHs heavier than 103M and a universal merger efficiency factor of pBH = 0.17 inferred from the IPTA data.

In the text
thumbnail Fig. A.1.

Distributions of the total signal strengths from 1, 10, 100, and 105 sources whose independent signal strengths have identical Ω 5 2 $ \Omega^{-\frac52} $ tails (left). Distribution of 𝒫circ for 1, 10 and 100 identical sources (right). The dashed line shows the distribution from a single source, the solid lines show the analytic estimates, and the distributions obtained from an explicit Monte Carlo analysis are shown by the points.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.