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/00046361/202346268  
Published online  02 August 2023 
Prospects for future binary black hole gravitational wave studies in light of PTA measurements
^{1}
Theoretical Physics Department, CERN, Geneva, Switzerland
^{2}
King’s College London, Strand, London WC2R 2LS, UK
^{3}
Keemilise ja Bioloogilise Füüsika Instituut, Rävala pst. 10, 10143 Tallinn, Estonia
email: juan.urrutia@kbfi.ee
^{4}
Departament of Cybernetics, Tallinn University of Technology, Akadeemia tee 21, 12618 Tallinn, Estonia
^{5}
Dipartimento di Fisica e Astronomia, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{6}
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
Received:
28
February
2023
Accepted:
3
June
2023
NANOGrav and other Pulsar Timing Arrays (PTAs) have discovered a commonspectrum process in the nHz range that may be due to gravitational waves (GWs): if so, they are likely to have been generated by black hole (BH) binaries with total masses > 10^{9} M_{⊙}. Using the Extended PressSchechter formalism to model the galactic halo mass function and a simple relation between the halo and BH masses suggests that these binaries have redshifts z = 𝒪(1) and mass ratios ≳10, and that the GW signal at frequencies above 𝒪(10) nHz may be dominated by relatively few binaries that could be distinguished experimentally and would yield observable circular polarization. Extrapolating the model to higher frequencies indicates that future GW detectors such as LISA and AEDGE could extend the PTA observations to lower BH masses ≳10^{3} M_{⊙}.
Key words: black hole physics / gravitational waves / cosmology: theory
© The Authors 2023
Open 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 𝒪(10^{6} − 10^{10}) 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 intermediatemass black holes (IMBHs) with masses in the range 𝒪(10^{3} − 10^{5}) 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 wellunderstood paths for forming compact stellarmass 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 commonspectrum 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 HellingsDowns 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 blackhole 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 commonspectrum process measured by the PTAs 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 nextgeneration 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 − 10^{2}) Hz and the projected AIONkm and AEDGE atom interferometers (Badurina et al. 2020; ElNeaj 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 > 𝒪(10^{9}) M_{⊙} that would probably be responsible for most of the PTA signal, LISA would observe a merger signal from binary black holes weighing 𝒪(10^{3} − 10^{9}) 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 𝒪(10^{3} − 10^{6}) 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, G_{N} = 1.
2. BH merger rate
As a first step towards estimating the BH merger rate, we use the Extended PressSchechter (EPS) formalism (Press & Schechter 1974; Bond et al. 1991) to calculate the galactic halo mass function,
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 M_{1} to merge with another one of mass, M_{2}, at some redshift, z, and become a halo of mass M_{f} = M_{1} + M_{2} (Lacey & Cole 1993):
Defining the halo merger rate kernel:
the halo merger rate is
As discussed by Benson et al. (2005), the merger rate kernel Q(M_{1}, M_{2}, t) in the EPS formalism is not symmetric in the exchange of M_{1} and M_{2}. However, it was shown by Erickcek et al. (2006) that taking M_{1} < M_{2} in the merger kernel aptly approximates the result of the BensonKamionkowskiHassani (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 p_{occ}(mM, z) that a BH of mass m occupies a halo of mass M, as well as the probability p_{merg}(m_{1}, m_{2}) that the galactic merger leads to a merger of their central BHs. The resulting BH merger rate can be expressed as
We assume the following simple redshiftdependent relation between the halo mass and the BH mass (Barkana & Loeb 2001; Wyithe & Loeb 2003):
where Δ_{c}(z) = 18π^{2} + 82[Ω_{M}(z)−1] − 39[Ω_{M}(z)−1]^{2} is the critical overdensity at virialization, which corresponds to p_{occ}(mM) = p_{occ}(m)δ(m − m_{BH}(M)). It is generally expected that SMBHs inhabit most large galaxies, and Xray observations (Miller et al. 2015) constrain the SMBH occupation fraction p_{occ}(m) to be > 20% for early galaxies with lower stellar masses 10^{7} < M_{*}/M_{⊙} < 10^{10}. 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 p_{occ} would shift this range towards lower redshifts when compared to a constant p_{occ}, increasing the observability of BH mergers.
In view of the paucity of information about the merger probability p_{merg}(m_{1}, m_{2}), we model p_{merg} by a constant for the sake of simplicity. After the galactic merger, various dynamical mechanisms must decrease the size of the SMBH binary below subparsec scales in order for a GW emissiondriven merger to take place (Begelman et al. 1980). The crossing of the final parsec is determined mostly by binary hardening via stellar losscone scattering (Merritt 2013), which can currently be considered the largest source of uncertainty in p_{merg}. The numerical study of Kelley et al. (2017a), based on a population of 10^{6} − 10^{10} M_{⊙} 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:
where p_{BH} ≡ p_{occ}(m_{1})p_{occ}(m_{2})p_{merg} 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 p_{BH} discussed above. We see that mergers with total masses M > 10^{9} M_{⊙} occur typically at z = 𝒪(1), those with total masses M ∈ (10^{6}, 10^{9}) M_{⊙} occur typically at z = 𝒪(3), and those with total masses M ∈ (10^{3}, 10^{6}) M_{⊙} occur typically at z = 𝒪(5).
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 p_{BH} discussed in the text. 
In this paper, we use the IPTA measurement (Antoniadis et al. 2022) to normalize p_{BH} for large BH masses, which circumvents the astrophysical uncertainties related to p_{occ} and p_{merg}, and extrapolate it to smaller BH masses. We also comment on the potential mass dependence of p_{BH}, which is relevant for the GW phenomenology when we extrapolate from the frequency band relevant for PTAs to those to be explored by nextgeneration 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):
where denotes the optimal amplitude of the Fourier transform of the GW strain, ,
is the differential merger rate and V_{c} denotes the comoving volume available at redshift z. We compute using the inspiralmergerringdown template (Ajith et al. 2008):
where the redshiftdependent chirp mass, ℳ_{z}, is given in terms of the total mass, M, and the symmetric mass ratio η ≡ m_{1}m_{2}/M^{2} of the binary by , and D_{L} denotes the luminosity distance of the binary. This implies the canonical scaling during the inspiral phase (Phinney 2001). The frequencies f_{merg}, f_{ring}, f_{cut}, and σ are parameterised as where a_{j}, b_{j}, and c_{j} 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 T_{tot}/T_{GW}, where T_{tot}, T_{GW} 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 (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 p_{BH} = 1, namely, the assumption that each galactic merger would produce a BH merger with m_{2} > m_{1} > 10^{3} M_{⊙}. This naive assumption leads to a spectrum exceeding that observed in the PTA band.
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 10^{3} M_{⊙} (p_{BH} = 1). The colored bands show the spectra from SMBHs heavier than 10^{9} M_{⊙} (dark blue), from SMBHs in the range (10^{6} M_{⊙}, 10^{9} M_{⊙}); light blue), and from IMBHs with masses in the range 10^{3} − 10^{6} M_{⊙} (green), assuming a universal efficiency factor . The shaded regions show the prospective sensitivities of IPTA30 (Chen et al. 2017; Kaiser & McWilliams 2021), LISA (LISA Collaboration 2017), AEDGE (ElNeaj et al. 2020; Badurina et al. 2021), AIONkm (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 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 p_{BH} so that the amplitude of the total GW spectrum from BH binaries matches the IPTA observations in the nHz range, which requires . 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 > 10^{9} M_{⊙} SMBH binaries, consistent with earlier studies (Wyithe & Loeb 2003; Sesana et al. 2004, 2008, 2009; Enoki et al. 2004; Kelley et al. 2017a; IzquierdoVillalba 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 ≲10^{10}, 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 ≳10^{4} 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 signaltonoise ratio of SNR = 8, which is similar to that expected for SKA (Janssen et al. 2015).
Fig. 3. Distributions of the binary parameters for a sampling of the BH binary population with 0 < z < 3, 1 nHz < f < 30nHz and . 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 ∈ (10^{6}, 10^{9}) M_{⊙}, shown in light blue, and the contribution in the AEDGE band is dominated by IMBHs with masses in the range of M ∈ (10^{3}, 10^{6}) M_{⊙}, shown in green. The lowfrequency GW emissions from the early stages of the inspirals of these binaries also contribute in the lowerfrequency bands, but these contributions are subdominant compared to the signals from highermass SMBH mergers^{4}.
We recall that the estimated spectrum in Fig. 2 assumes a constant p_{BH}. Clearly, the relative contribution from binaries with lighter BH masses < 10^{9} M_{⊙} could be enhanced (reduced) by increasing (decreasing) p_{BH} for these binaries. In order for mergers of BHs with masses in the range (10^{6}, 10^{9}) M_{⊙} to dominate the PTA signal p_{BH} ≃ 1 would be required in this mass range, but even this maximal enhancement would be insufficient for mergers with M < 10^{6} M_{⊙} 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 p_{BH} ≃ 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 p_{BH} 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)
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 is^{5}
We show in Fig. 4 the expected numbers of binaries with total masses of M ≥ M_{min} 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 10^{5} (10^{2}) binaries with M > 3 × 10^{9} M_{⊙} 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 higherfrequency bins, may be comprised of a limited number of nearly monochromatic signals from heavy SMBH binaries ( > 10^{9} M_{⊙}), implying sizeable fluctuations around the smooth spectrum that would be obtained in the limit of a large population of heavy inspiralling binaries.
Fig. 4. Expected number of SMBH binaries heavier than M_{min} and with symmetric mass ratio η > 0.01 emitting GWs in the indicated frequency range for p_{BH} = 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 (f_{j}, f_{j + 1}) populated with N(f_{j}) binaries. A realization of the GW background from each bin can then be obtained as (see Eq. (A.6)):
where the contribution from an individual binary emitting in this frequency band is given by
and θ ≡ {ℳ_{z},zη,f} are parameters describing the binary. The number of binaries N(f_{j}) in this frequency bin is drawn from a Poisson distribution with the expectation value and the binary parameters are generated randomly according to the distribution
This is the expression that was used to calculate the distributions of binary parameters shown in Fig. 3 for a sampling with p_{BH} = 0.17 of 4 × 10^{5} 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 , so that 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 emitting at a frequency f, which can be expressed via two oneparameter 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 expectation in higher frequency bins, while a few bins display sharp peaks.
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 , 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, p_{BH} = 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 orderofmagnitude estimate, we focus on the closest binaries and ignore the redshift dependence. In this case, , by Eqs. (8) and (10), and the probability of finding an event at D_{L} is . Thus , where, in this simplified approach, D_{L,max} is some large luminosity distance at which gets suppressed, and D_{L,min} is the distance to the nearest possible source. On the other hand, . Therefore, we expect the mean of the GW signal to be determined by faraway sources, while the variance is set by a few closeby binaries. Moreover, we can estimate from P^{(1)}(D_{L}) that has a relatively flat powerlaw tail at large values,
Since the closest distance to massive BH binaries is constrained, this tail will be cut off at . 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 (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., ), the distribution of Ω_{GW} does not converge to a Gaussian because it retains its 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 widthtomean ratio scales as:
when , as is shown in Appendix A. This scaling is slower than the typical 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; SatoPolito & 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
where
is the maximal GW strain from an inspiralling binary, and the gravitational Stokes parameters are defined by
The fractional amount of circular polarization of the SGWB can be characterized by the quantity
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 (N_{dom} ≳ 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)
The lowerleft 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 > 10^{9} M_{⊙} 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 lowerright 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 p_{BH} is not universal. By suppressing p_{BH}, for instance, at M > 10^{9} M_{⊙} and enhancing p_{BH} for 10^{9} M_{⊙} > M > 10^{6} M_{⊙}, we can accommodate the IPTA commonspectrum effect. In this case, the GW energy spectrum Ω_{GW} becomes smoother and approaches the naive behavior while suppressing the circular polarization 𝒫_{circ}. Moreover, in this case, p_{BH} 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 slowlyevolving 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
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 S_{n}(f) as
The optimal timedependent 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
where
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 p_{det} 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 S_{n} include the foregrounds from stellar mass BH binaries and white dwarf binaries (see e.g., Lewicki & Vaskonen 2021). Furthermore, we used SNR_{c} = 8 for the detection threshold and 𝒯 = 1 yr for the observation time.
In the monochromatic limit, Eq. (24) is simplified to . Using this, we can convert the prospected noise S_{n}(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 nearmonochromatic 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 nearmonochromatic sources because the coalescence time of binaries heavier than 10^{3} M_{⊙} is less than a day when they enter the AEDGE sensitivity window^{6}. 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} = 10^{4} M_{⊙} reflecting the binary population, that increases towards lower masses, and the cutoff at 10^{3} M_{⊙}. 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.
Fig. 6. Expected numbers of nearmonochromatic 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 10^{3} M_{⊙} and a universal merger efficiency factor of p_{BH} = 0.17 inferred from IPTA data, as discussed in the text. 
We find that, at z < 10, there are a total of 𝒪(10^{7}) nearmonochromatic sources in the frequency range of 10^{−5} Hz < f < 0.1 Hz^{7}. 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 p_{BH} 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} > 10^{4} M_{⊙}, while AEDGE can see lighter binaries, ℳ_{0} < 2 × 10^{3} M_{⊙}. LISA could (in principle) detect even heavier mergers, namely, of ℳ_{0} ≲ 10^{8} M_{⊙}, but since such mergers are so rare it is unlikely that LISA will see any mergers above 10^{6} M_{⊙}. 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} ∼ 10^{4} 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 AIONkm and ET we find that the total expected number of events per year is less than three.
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 10^{3} M_{⊙} and a universal merger efficiency factor of p_{BH} = 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 multimessenger signals.
We caution that there are observational uncertainties in the low mass cutoff 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 lowmass region and are compatible with the assumed cutoff, although a heavier mass cutoff 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 cutoff 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 cutoff 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 commonprocess signal based on a simulation of massive BH mergers. The magnitude of the signal depends on the merger probability, p_{BH}, 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 massindependent value of p_{BH} to the PTA signal, finding p_{BH} = 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 > 10^{9} M_{⊙} and about 10% from masses < 10^{9} M_{⊙}. 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 massindependent value of p_{BH} as for the PTA signal, there would be observable signals from mergers with total masses ∈(10^{3}, 10^{6}) M_{⊙} in the LISA experiment and from mergers with total masses ∈(10^{3}, 10^{5}) M_{⊙} in the AEDGE experiment. Data from these experiments will be able to check the accuracy of the constancy of p_{BH}. In principle, their signals could be even larger than our estimates if p_{BH} 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 p_{BH} and, hence, a lower event rate cannot be excluded.
However, alternative scenarios for quasar optical variability may also be considered, such as intrinsic variability in the accretion disc (King et al. 2013).
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., BlancoPillado et al. (2021).
We note in passing that binaries capable of generating the PTA nHz background could not explain the yearlike quasiperiodicity seen in blazars (Holgado et al. 2018).
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 p_{BH} = 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 powerlaw 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.
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.
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łodowskaCurie grant agreement No. 101065736.
References
 Ajith, P., Babak, S., Chen, Y., et al. 2008, Phys. Rev. D, 77, 104017 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadis, J., Arzoumanian, Z., Babak, S., et al. 2022, MNRAS, 510, 4873 [NASA ADS] [CrossRef] [Google Scholar]
 Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Arzoumanian, Z., Baker, P. T., Blecha, L., et al. 2023, ApJL, submitted [arXiv:2301.03608] [Google Scholar]
 Ashoorioon, A., Rezazadeh, K., & Rostami, A. 2022, Phys. Lett. B, 835, 137542 [NASA ADS] [CrossRef] [Google Scholar]
 Badurina, L., Bentine, E., Blas, D., et al. 2020, JCAP, 05, 011 [CrossRef] [Google Scholar]
 Badurina, L., Buchmueller, O., Ellis, J., et al. 2021, Phil. Trans. A. Math. Phys. Eng. Sci., 380, 20210060 [Google Scholar]
 Barkana, R., & Loeb, A. 2001, Phys. Rept., 349, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022, ApJ, 941, 119 [CrossRef] [Google Scholar]
 Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
 Belgacem, E., & Kamionkowski, M. 2020, Phys. Rev. D, 102, 023004 [NASA ADS] [CrossRef] [Google Scholar]
 Benson, A. J., Kamionkowski, M., & Hassani, S. H. 2005, MNRAS, 357, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Berti, E., Cardoso, V., & Will, C. M. 2006, Phys. Rev. D, 73, 064030 [NASA ADS] [CrossRef] [Google Scholar]
 BlancoPillado, J. J., Olum, K. D., & Wachter, J. M. 2021, Phys. Rev. D, 103, 103512 [NASA ADS] [CrossRef] [Google Scholar]
 Blas, D., & Jenkins, A. C. 2022, Phys. Rev. Lett., 128, 101103 [NASA ADS] [CrossRef] [Google Scholar]
 Blasi, S., Brdar, V., & Schmitz, K. 2021, Phys. Rev. Lett., 126, 041305 [NASA ADS] [CrossRef] [Google Scholar]
 Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53 [Google Scholar]
 Buchmuller, W., Domcke, V., & Schmitz, K. 2020, Phys. Lett. B, 811, 135914 [NASA ADS] [CrossRef] [Google Scholar]
 Chadayammuri, U., Bogdan, A., Ricarte, A., & Natarajan, P. 2023, ApJ, 946, 51 [CrossRef] [Google Scholar]
 Chen, S., Middleton, H., Sesana, A., Del Pozzo, W., & Vecchio, A. 2017, MNRAS, 468, 404 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970 [NASA ADS] [CrossRef] [Google Scholar]
 Conneely, C., Jaffe, A. H., & Mingarelli, C. M. F. 2019, MNRAS, 487, 562 [NASA ADS] [CrossRef] [Google Scholar]
 De Luca, V., Franciolini, G., & Riotto, A. 2021, Phys. Rev. Lett., 126, 041303 [NASA ADS] [CrossRef] [Google Scholar]
 De Paolis, F., Ingrosso, G., & Nucita, A. A. 2002, A&A, 388, 470 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Paolis, F., Ingrosso, G., & Nucita, A. A. 2004, A&A, 426, 379 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dodelson, S. 2003, Modern Cosmology (Amsterdam: Academic Press) [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [Google Scholar]
 ElNeaj, Y. A., Alpigiani, C., AmairiPyka, S., et al. 2020, EPJ Quant. Technol., 7, 6 [CrossRef] [Google Scholar]
 Ellis, J., & Lewicki, M. 2021, Phys. Rev. Lett., 126, 041304 [NASA ADS] [CrossRef] [Google Scholar]
 Enoki, M., & Nagashima, M. 2007, Prog. Theor. Phys., 117, 241 [CrossRef] [Google Scholar]
 Enoki, M., Inoue, K. T., Nagashima, M., & Sugiyama, N. 2004, ApJ, 615, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Erickcek, A. L., Kamionkowski, M., & Benson, A. J. 2006, MNRAS, 371, 1992 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Finn, L. S., & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198 [NASA ADS] [CrossRef] [Google Scholar]
 Fragione, G., & Kocsis, B. 2018, Phys. Rev. Lett., 121, 161103 [Google Scholar]
 Fragione, G., & Loeb, A. 2023, ApJ, 944, 81 [CrossRef] [Google Scholar]
 Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Gerosa, D., Ma, S., Wong, K. W. K., et al. 2019, Phys. Rev. D, 99, 103004 [NASA ADS] [CrossRef] [Google Scholar]
 Goncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021, ApJ, 917, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74 [Google Scholar]
 Graham, P. W., Hogan, J. M., Kasevich, M. A., Rajendran, S., & Romani, R. W. 2017, ArXiv eprints [arXiv:1711.02225] [Google Scholar]
 Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
 Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Hogan, J. M., Johnson, D. M. S., Dickerson, S., et al. 2011, Gen. Rel. Grav., 43, 1953 [NASA ADS] [CrossRef] [Google Scholar]
 Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, MNRAS, 481, L74 [NASA ADS] [CrossRef] [Google Scholar]
 Hotinli, S. C., Kamionkowski, M., & Jaffe, A. H. 2019, Open J. Astrophys., 2, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Hütsi, G., Raidal, M., Urrutia, J., Vaskonen, V., & Veermäe, H. 2023, Phys. Rev. D, 107, 043502 [CrossRef] [Google Scholar]
 Iguchi, S., Okuda, T., & Sudou, H. 2010, ApJ, 724, L166 [NASA ADS] [CrossRef] [Google Scholar]
 IzquierdoVillalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2021, MNRAS, 509, 3488 [CrossRef] [Google Scholar]
 Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, PoS, AASKA14, 037 [Google Scholar]
 Kaiser, A. R., & McWilliams, S. T. 2021, CQG, 38, 055009 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, R., & Soda, J. 2016, Phys. Rev. D, 93, 062003 [NASA ADS] [CrossRef] [Google Scholar]
 Kelley, L. Z., Blecha, L., & Hernquist, L. 2017a, MNRAS, 464, 3131 [NASA ADS] [CrossRef] [Google Scholar]
 Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017b, MNRAS, 471, 4508 [NASA ADS] [CrossRef] [Google Scholar]
 King, O. G., Hovatta, T., MaxMoerbeck, W., et al. 2013, MNRAS, 436, 114 [Google Scholar]
 Kocsis, B., & Sesana, A. 2011, MNRAS, 411, 1467 [NASA ADS] [CrossRef] [Google Scholar]
 Kohri, K., & Terada, T. 2021, Phys. Lett. B, 813, 136040 [NASA ADS] [CrossRef] [Google Scholar]
 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]
 Lacey, C. G., & Cole, S. 1993, MNRAS, 262, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Lewicki, M., & Vaskonen, V. 2021, Eur. Phys. J. C, 83, 168 [Google Scholar]
 LIGO Scientific & Virgo Collaborations (Abbott, B. P., et al.) 2016, Phys. Rev. Lett., 116, 061102 [CrossRef] [PubMed] [Google Scholar]
 LIGO Scientific, Virgo & KAGRA Collaborations (Abbott, R., et al.) 2021, ArXiv eprints [arXiv:2111.03634] [Google Scholar]
 Lippai, Z., Frei, Z., & Haiman, Z. 2009, ApJ, 701, 360 [CrossRef] [Google Scholar]
 LISA Collaboration (AmaroSeoane, P., et al.) 2017, ArXiv eprints [arXiv:1702.00786] [Google Scholar]
 Merritt, D. 2013, CQG, 30, 244005 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98 [NASA ADS] [CrossRef] [Google Scholar]
 O’Neill, S., Kiehlmann, S., Readhead, A. C. S., et al. 2022, ApJ, 926, L35 [CrossRef] [Google Scholar]
 Peters, P. C. 1964, Phys. Rev., 136, B1224 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, P. C., & Mathews, J. 1963, Phys. Rev., 131, 435 [Google Scholar]
 Phinney, E. S. 2001, ArXiv eprints [arXiv:astroph/0108028] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
 Rieger, F. M., & Mannheim, K. 2000, A&A, 359, 948 [NASA ADS] [Google Scholar]
 Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, CQG, 29, 124013 [NASA ADS] [CrossRef] [Google Scholar]
 SatoPolito, G., & Kamionkowski, M. 2022, Phys. Rev. D, 106, 023004 [NASA ADS] [CrossRef] [Google Scholar]
 Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, ApJ, 611, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Sesana, A., Vecchio, A., & Volonteri, M. 2009, MNRAS, 394, 2255 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, S. R., Huerta, E. A., Gair, J. R., & McWilliams, S. T. 2016, ApJ, 817, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Valbusa Dall’Armi, L., Nishizawa, A., Ricciardone, A., & Matarrese, S. 2023, Phys. Rev. Lett., 131, 041401 [CrossRef] [Google Scholar]
 Valtonen, M. J., et al. 2008, Nature, 452, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Vaskonen, V., & Veermäe, H. 2021, Phys. Rev. Lett., 126, 051303 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [Google Scholar]
 Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nat. Rev. Phys., 3, 732 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y., Pardo, K., Chang, T.C., & Doré, O. 2021, Phys. Rev. D, 103, 084007 [NASA ADS] [CrossRef] [Google Scholar]
 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:
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
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 are independent and identically distributed. Thus, the statistical properties of Ω_{GW}(f) can be inferred from the distribution of , which is given by
where we define the following function:
and the spectral source density:
We see that, for circular inspiralling binaries, the fourparameter 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 (f_{j}, f_{j + 1}) can be obtained from
where the number of binaries N(f_{j}) is drawn from a Poisson distribution with the expected value,
and the are drawn from a random distribution with f ∈ (f_{j}, f_{j + 1}). In the following, we suppress f_{j}.
The moment generating function of Ω_{GW} is given by:
Thus, the cumulant generating function of Ω_{GW}(f_{j}), that is,
is proportional to the generating function of , implying that the nth cumulant of Ω_{GW} is
Therefore, the mean κ_{1}(Ω_{GW})≡⟨Ω_{GW}⟩ matches the expectation value given in Eq. (8), while the variance,
is divergent due to the long 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 .
The momentgenerating 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}(f_{j}) can be obtained by an inverse Laplace transform of the moment generating function given in Eq. (A.8). It can be expressed as
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 . 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):
where C is a constant. The second moment ⟨Ω^{2}⟩ diverges. We assume that the small Ω behavior is such that the first moment is finite. The 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):
In summary, the N → ∞ distribution asymptotes to
where denotes the expectation value of Ω and
is a universal function that is independent of the single event distribution – all information about the latter enters through and C. This limiting case reveals a few crucial features:

Since 𝒫(x)∼x^{−5/2} when x → ∞, the powerlaw behavior of the tail is preserved in the large limit:

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

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
So, the relative width will approach zero, but at a slower pace than when the central limit theorem applies.

Using the method of steepest decent, we find that
implying values smaller than the mean, that is , 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 longtailed distribution
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 if the signal consists of 10, 100, and 10^{5} sources. The mean is given by . 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 = 10^{5}, but predicts a slightly too wide peak for N ≤ 10^{5}. 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.
Fig. A.1. Distributions of the total signal strengths from 1, 10, 100, and 10^{5} sources whose independent signal strengths have identical 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
which depends on the distribution and number N of sources. We denote the contributions from individual sources by V_{i}, I_{i} given by Eqs. (1820). 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
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
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}h_{i}^{2}, but not over signal strengths, so that
where p_{i} ≡ (4/5)h_{i}^{2}/⟨I⟩_{θ} ∈ [0, 1] is the fractional contribution to the signal from source i and 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 p_{i}, which vary only in the range [0, 1]. We find:
We note that N⟨p_{i}^{2}⟩ takes values in the range [1/N, 1] and could be interpreted as a measure of the effective number of dominant sources 1/N_{dom}. For instance, in the limiting case where the signal is sourced by N binaries that contribute equally, that is, p_{i} = 1/N, we would obtain
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),
The corrections arising from higher powers of δI/⟨I⟩_{θ} are suppressed by increasing powers of 1/N.
All Figures
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 p_{BH} discussed in the text. 

In the text 
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 10^{3} M_{⊙} (p_{BH} = 1). The colored bands show the spectra from SMBHs heavier than 10^{9} M_{⊙} (dark blue), from SMBHs in the range (10^{6} M_{⊙}, 10^{9} M_{⊙}); light blue), and from IMBHs with masses in the range 10^{3} − 10^{6} M_{⊙} (green), assuming a universal efficiency factor . The shaded regions show the prospective sensitivities of IPTA30 (Chen et al. 2017; Kaiser & McWilliams 2021), LISA (LISA Collaboration 2017), AEDGE (ElNeaj et al. 2020; Badurina et al. 2021), AIONkm (Badurina et al. 2020, 2021) and ET (Sathyaprakash et al. 2012). 

In the text 
Fig. 3. Distributions of the binary parameters for a sampling of the BH binary population with 0 < z < 3, 1 nHz < f < 30nHz and . 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 
Fig. 4. Expected number of SMBH binaries heavier than M_{min} and with symmetric mass ratio η > 0.01 emitting GWs in the indicated frequency range for p_{BH} = 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 
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 , 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, p_{BH} = 0.17 is assumed. 

In the text 
Fig. 6. Expected numbers of nearmonochromatic 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 10^{3} M_{⊙} and a universal merger efficiency factor of p_{BH} = 0.17 inferred from IPTA data, as discussed in the text. 

In the text 
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 10^{3} M_{⊙} and a universal merger efficiency factor of p_{BH} = 0.17 inferred from the IPTA data. 

In the text 
Fig. A.1. Distributions of the total signal strengths from 1, 10, 100, and 10^{5} sources whose independent signal strengths have identical 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.