11institutetext:  Theoretical Physics Department, CERN, Geneva, Switzerland
22institutetext: King’s College London, Strand, London, WC2R 2LS, United
Kingdom 33institutetext: Keemilise ja Bioloogilise Füüsika Instituut, Rävala
pst. 10, 10143 Tallinn, Estonia 44institutetext: Departament of Cybernetics,
Tallinn University of Technology, Akadeemia tee 21, 12618 Tallinn, Estonia
55institutetext: Dipartimento di Fisica e Astronomia, Università degli Studi
di Padova, Via Marzolo 8, 35131 Padova, Italy 66institutetext: Istituto
Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova,
Italy

# Prospects for future binary black hole gravitational wave studies in light
of PTA measurements

John Ellis , 112233 Malcolm Fairbairn , 22 Gert Hütsi , 33 Martti Raidal , 33
Juan Urrutia<EMAIL_ADDRESS>
Ville Vaskonen , 335566 Hardi Veermäe , 33

(KCL-PH-TH/2023-04, CERN-TH-2023-008, AION-REPORT/2023-1)

NANOGrav and other Pulsar Timing Arrays (PTAs) have discovered a common-
spectrum 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_{\odot}$. Using the Extended Press-
Schechter 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={\cal O}(1)$ and mass ratios $\gtrsim 10$, and that the GW signal
at frequencies above ${\cal O}(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 $\gtrsim 10^{3}M_{\odot}$.

###### Key Words.:

Gravitational waves – quasars: supermassive black holes – Stars: black holes

##  1 Introduction

The discovery (Abbott et al. 2016) of gravitational waves (GWs) by the LIGO
and Virgo experiments, in the ${\cal O}(100)$ Hz range of frequencies
generated by mergers of black holes (BHs) with masses of ${\cal
O}(10-100)M_{\odot}$ (Abbott et al. 2021), has opened a new window onto both
astrophysics and cosmology. Supermassive black holes (SMBHs) with masses of
${\cal O}(10^{6}-10^{10})M_{\odot}$ 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 (Akiyama et al. 2019, 2022). However,
information on intermediate-mass black holes (IMBHs) with masses in the range
${\cal O}(10^{3}-10^{5})M_{\odot}$ (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 et al. 2008; Boroson &
Lauer 2009; Iguchi et al. 2010; Graham et al. 2015; O’Neill et al. 2022;
Kovačević et al. 2022) 111However, alternative scenarios for quasar optical
variability may also be considered, such as intrinsic variability in the
accretion disc (King et al. 2013).. 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), 222A 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). but they have not (yet)
detected the Hellings-Downs quadrupolar signature that is characteristic of
GWs (Hellings & Downs 1983) 333We 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)..

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^{+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
$\alpha=-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
${\cal O}(10^{-4}-10^{-2})$ Hz, the proposed ET laser interferometer
(Sathyaprakash et al. 2012), whose peak sensitivity is at frequencies ${\cal
O}(1-10^{2})$ 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 ${\cal O}(10^{-2}-1)$ Hz. We show that, while none of these
detectors would observe a signal from the binary black holes weighing $>{\cal
O}(10^{9})M_{\odot}$ that would probably be responsible for most of the PTA
signal, LISA would observe a merger signal from binary black holes weighing
${\cal O}(10^{3}-10^{9})M_{\odot}$ 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 ${\cal O}(10^{3}-10^{6})M_{\odot}$ 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
Press-Schechter (EPS) formalism (Press & Schechter 1974; Bond et al. 1991) to
calculate the galactic halo mass function,

$\frac{{\rm d}n(M,t)}{{\rm d}\ln
M}=\frac{\rho_{0}}{M}\sqrt{\frac{2}{\pi}}\frac{{\rm d}\ln\sigma}{{\rm d}\ln
M}\frac{\delta_{c}(z)}{\sigma(M)}e^{-\frac{\delta_{c}(z)^{2}}{2\sigma^{2}(M)}}\,,$
(1)

where $\rho_{0}$ is the background matter density today, $\sigma^{2}(M)$ is
the variance of the matter fluctuations and $\delta_{c}(z)$ is the critical
overdensity for collapse. The latter is given by $\delta_{c}(z)\approx
1.686/D(z)$ where $D(z)$ is the linear growth function (Dodelson 2003). We
calculate $\sigma^{2}(M)$ from the cold dark matter power spectrum with the
Planck 2018 cosmological parameters (Aghanim et al. 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):

$\displaystyle\frac{{\rm d}p(M_{1},M_{2},t)}{{\rm d}t{\rm
d}M_{2}}=\frac{1}{M_{f}}\sqrt{\frac{2}{\pi}}\left|\frac{\dot{\delta}_{c}}{\delta_{c}}\right|\frac{{\rm
d}\ln{\sigma}}{{\rm
d}\ln{M_{f}}}\left[1-\frac{\sigma^{2}(M_{f})}{\sigma^{2}(M_{1})}\right]^{-\frac{3}{2}}$
(2)
$\displaystyle\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]\,.$

Defining the halo merger rate kernel:

$Q(M_{1},M_{2},t)\equiv\frac{{\rm d}p(M_{1},M_{2},t)}{{\rm d}t{\rm
d}M_{2}}\left[\frac{{\rm d}n(M_{2},t)}{{\rm d}M_{2}}\right]^{-1},$ (3)

the halo merger rate is

$\frac{{\rm d}R_{h}}{{\rm d}M_{1}{\rm d}M_{2}}=\frac{{\rm d}n(M_{1},t)}{{\rm
d}M_{1}}\frac{{\rm d}n(M_{2},t)}{{\rm d}M_{2}}Q(M_{1},M_{2},t)\,.$ (4)

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
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.

Figure 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_{\rm BH}$ discussed in the text.

In order to estimate the merger rate of central BHs, we must estimate the
probability $p_{\rm occ}(m|M,z)$ that a BH of mass $m$ occupies a halo of mass
$M$, as well as the probability $p_{\rm 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

$\displaystyle\frac{{\rm d}R_{\rm BH}}{{\rm d}m_{1}{\rm d}m_{2}}=$
$\displaystyle\int{\rm d}M_{1}{\rm d}M_{2}\,p_{\rm merg}(m_{1},m_{2})$ (5)
$\displaystyle\times p_{\rm occ}(m_{1}|M_{1},z)p_{\rm
occ}(m_{2}|M_{2},z)\frac{{\rm d}R_{h}}{{\rm d}M_{1}{\rm d}M_{2}}\,.$

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

$\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}}\,,$ (6)

where $\Delta_{c}(z)=18\pi^{2}+82[\Omega_{M}(z)-1]-39[\Omega_{M}(z)-1]^{2}$ is
the critical overdensity at virialization, which corresponds to $p_{\rm
occ}(m|M)=p_{\rm occ}(m)\delta(m-m_{\rm BH}(M))$. It is generally expected
that SMBHs inhabit most large galaxies, and X-ray observations (Miller et al.
2015) constrain the SMBH occupation fraction $p_{\rm occ}(m)$ to be $>20\%$
for early galaxies with lower stellar masses $10^{7}<M_{*}/M_{\odot}<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\lesssim z\lesssim 3$. The estimated $z$ dependence of $p_{\rm occ}$
would shift this range towards lower redshifts when compared to a constant
$p_{\rm occ}$, increasing the observability of BH mergers.

Figure 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_{\odot}$ ($p_{\rm BH}=1$). The colored bands show the spectra
from SMBHs heavier than $10^{9}M_{\odot}$ (dark blue), from SMBHs in the range
$(10^{6}M_{\odot},10^{9}M_{\odot})$ (light blue), and from IMBHs with masses
in the range $10^{3}-10^{6}M_{\odot}$ (green), assuming a universal efficiency
factor $p_{\rm BH}=0.17^{+0.18}_{-0.08}$. The shaded regions show the
prospective sensitivities of IPTA30 (Chen et al. 2017; Kaiser & McWilliams
2021), LISA (Amaro-Seoane et al. 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 view of the paucity of information about the merger probability $p_{\rm
merg}(m_{1},m_{2})$, we model $p_{\rm 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 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 $p_{\rm merg}$. The numerical study of Kelley et al.
(2017a), based on a population of $10^{6}-10^{10}M_{\odot}$ 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:

$\frac{{\rm d}R_{\rm BH}}{{\rm d}m_{1}{\rm d}m_{2}}\approx p_{\rm
BH}\frac{{\rm d}M_{1}}{{\rm d}m_{1}}\frac{{\rm d}M_{2}}{{\rm
d}m_{2}}\frac{{\rm d}R_{h}}{{\rm d}M_{1}{\rm d}M_{2}}\,,$ (7)

where $p_{\rm BH}\equiv p_{\rm occ}(m_{1})p_{\rm occ}(m_{2})p_{\rm 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_{\rm BH}$ discussed
above. We see that mergers with total masses $M>10^{9}M_{\odot}$ occur
typically at $z={\cal O}(1)$, those with total masses
$M\in(10^{6},10^{9})M_{\odot}$ occur typically at $z={\cal O}(3)$, and those
with total masses $M\in(10^{3},10^{6})M_{\odot}$ occur typically at $z={\cal
O}(5)$.

In this paper, we use the IPTA measurement (Antoniadis et al. 2022) to
normalize $p_{\rm BH}$ for large BH masses, which circumvents the
astrophysical uncertainties related to $p_{\rm occ}$ and $p_{\rm merg}$, and
extrapolate it to smaller BH masses. We also comment on the potential mass
dependence of $p_{\rm BH}$, 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

Figure 3: Distributions of the binary parameters for a sampling of the BH
binary population with $0<z<3$, $1\,{\rm nHz}<f<30{\rm nHz}$ and
$\mathcal{M}_{0}>10^{9}(f/{\rm nHz})^{-\frac{4}{3}}$. The symmetric mass ratio
is integrated with the lower bound $\eta>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$.

###  3.1 The GW energy spectrum

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

$\Omega_{\rm GW}(f)\equiv\frac{1}{\rho_{\rm c}}\frac{{\rm d}\rho_{\rm
GW}}{{\rm d}\ln
f}=\int\text{d}\lambda\,\frac{2\pi}{5}\frac{f^{3}|\tilde{h}(f)|^{2}}{\rho_{c}}\,,$
(8)

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

${\rm d}\lambda={\rm d}m_{1}{\rm d}m_{2}\,\frac{{\rm d}z}{1+z}\frac{{\rm
d}V_{c}}{{\rm d}z}\frac{{\rm d}R_{\rm BH}(t)}{{\rm d}m_{1}{\rm d}m_{2}}\,,$
(9)

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

$\displaystyle|\tilde{h}(f)|$
$\displaystyle=\sqrt{\frac{5}{24}}\frac{\mathcal{M}_{z}^{\frac{5}{6}}}{\pi^{\frac{2}{3}}D_{L}}$
(10) $\displaystyle\times\begin{cases}f^{-\frac{7}{6}}&f<f_{\rm merg}\\\
f_{\rm merg}^{-\frac{1}{2}}f^{-\frac{2}{3}}&f_{\rm merg}\leq 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}\leq f<f_{\rm cut}\,,\end{cases}$

where the redshift-dependent chirp mass, $\mathcal{M}_{z}$, is given in terms
of the total mass , $M$, and the symmetric mass ratio $\eta\equiv
m_{1}m_{2}/M^{2}$ of the binary by
$\mathcal{M}_{z}\equiv(1+z)\mathcal{M}_{0}=(1+z)M\eta^{\frac{3}{5}}$, and
$D_{L}$ denotes the luminosity distance of the binary. This implies the
canonical $\Omega_{\rm GW}\propto f^{\frac{2}{3}}$ scaling during the inspiral
phase (Phinney 2001). The frequencies $f_{\rm merg}$, $f_{\rm ring}$, $f_{\rm
cut}$, and $\sigma$ are parameterised as
$f_{j}=\eta^{\frac{5}{3}}(a_{j}\eta^{2}+b_{j}\eta+c_{j})/(\pi\mathcal{M}_{z})$
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_{\rm tot}/T_{\rm GW}$, where
$T_{\rm tot}$, $T_{\rm 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
$(1+73/24e^{2}+37/96e^{4})/(1-e^{2})^{\frac{7}{2}}$ (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_{\rm BH}=1$, namely, the
assumption that each galactic merger would produce a BH merger with
$m_{2}>m_{1}>10^{3}M_{\odot}$. This naive assumption leads to a spectrum
exceeding that observed in the PTA band.

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 $\Omega_{\rm GW}\propto
f^{\frac{2}{3}}$ tail at low frequencies arises from nearly monochromatic
signals generated by inspiralling BH binaries. Because of this, $\Omega_{\rm
GW}$ will fluctuate depending on the specific realization of the binary
population, so the spectrum in Fig. 2 should be interpreted as the mean
$\Omega_{\rm GW}$ obtained by averaging over many potential realizations of
the binary population. We return to the probability distribution of
$\Omega_{\rm GW}$ in the next subsection and in Appendix A.

The grey band in Fig. 2 was obtained by choosing a fixed value of $p_{\rm BH}$
so that the amplitude of the total GW spectrum from BH binaries matches the
IPTA observations in the nHz range, which requires $p_{\rm
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>10^{9}M_{\odot}$ SMBH binaries, consistent with earlier studies (Wyithe &
Loeb 2003; Sesana et al. 2004; Enoki et al. 2004; Sesana et al. 2008, 2009;
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 $\lesssim 10^{10}$, and the redshifts are
typically ${\cal O}(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 $\gtrsim 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 signal-to-noise ratio of SNR = 8, which is similar to that expected
for SKA (Janssen et al. 2015).

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\in(10^{6},10^{9})M_{\odot}$, shown in light blue, and the contribution in
the AEDGE band is dominated by IMBHs with masses in the range of
$M\in(10^{3},10^{6})M_{\odot}$, 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 mergers. 444In 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_{\rm BH}=1$, but a possible
improvement of the Roman sensitivity might enable the detection of SMBH binary
signals at frequencies $\in(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\sim 10^{-6}$ Hz is well below that
(see Sec. 3.2) and consists of a few nearly monochromatic sources.

We recall that the estimated spectrum in Fig. 2 assumes a constant $p_{\rm
BH}$. Clearly, the relative contribution from binaries with lighter BH masses
$<10^{9}M_{\odot}$ could be enhanced (reduced) by increasing (decreasing)
$p_{\rm BH}$ for these binaries. In order for mergers of BHs with masses in
the range $(10^{6},10^{9})M_{\odot}$ to dominate the PTA signal $p_{\rm
BH}\simeq 1$ would be required in this mass range, but even this maximal
enhancement would be insufficient for mergers with $M<10^{6}M_{\odot}$ 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_{\rm BH}\simeq 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_{\rm 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)

$\frac{{\rm d}N}{{\rm d}\mathcal{M}_{z}{\rm d}\eta{\rm d}z{\rm d}\ln
f}=\frac{{\rm d}\lambda}{{\rm d}\mathcal{M}_{z}{\rm d}\eta{\rm
d}z}\left|\frac{{\rm d}\tau}{{\rm d}\ln f}\right|\,,$ (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$ is 555We 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).

$\tau(f)=\frac{5}{256\pi^{\frac{8}{3}}\mathcal{M}_{z}^{\frac{5}{3}}}f^{-\frac{8}{3}}\,.$
(12)

Figure 4: Expected number of SMBH binaries heavier than $M_{\rm min}$ and with
symmetric mass ratio $\eta>0.01$ emitting GWs in the indicated frequency range
for $p_{\rm 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.

We show in Fig. 4 the expected numbers of binaries with total masses of $M\geq
M_{\rm min}$ and mass ratios $q\geq 0.01$ that emit GWs in several frequency
bands in the range $1\,{\rm nHz}<f<30\,{\rm 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\times 10^{9}M_{\odot}$
generate 80% of the calculated GW signal in the $1\,{\rm nHz}<f<3\,{\rm nHz}$
($10\,{\rm nHz}<f<30\,{\rm nHz}$) frequency band, whereas 50% of the
calculated GW signal in the $3\,{\rm nHz}<f<10\,{\rm 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
$(>10^{9}M_{\odot})$, implying sizeable fluctuations around the smooth
$f^{\frac{2}{3}}$ spectrum that would be obtained in the limit of a large
population of heavy inspiralling binaries.

Figure 5: GW energy spectrum $\Omega_{\rm GW}$ and the fractional circular
polarization $\mathcal{P}_{\rm 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 $\Omega_{\rm GW}$ and
$\mathcal{P}_{\rm 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\propto\Omega^{-\frac{5}{2}}$, and in the lower
right panel, the dashed curves show Gaussian fits to these distributions of
width $\sigma=0.13$ (blue) and $\sigma=0.23$ (orange). In this figure, $p_{\rm
BH}=0.17$ is assumed.

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. (32)):

$\Omega_{\rm
GW}(f_{j})=\frac{1}{\ln(f_{j+1}/f_{j})}\sum_{k=1}^{N(f_{j})}\Omega_{\rm
GW}^{(1)}(\@vec{\theta}_{k})\,,$ (13)

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

$\Omega_{\rm GW}^{(1)}(\@vec{\theta})=\left|\frac{{\rm d}\tau}{{\rm d}\ln
f}\right|^{-1}\frac{2\pi}{5}\frac{f^{3}|\tilde{h}(f)|^{2}}{\rho_{c}}\,,$ (14)

and $\@vec{\theta}\equiv\\{\mathcal{M}_{z},z,\eta,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
$\bar{N}(f_{j})\equiv\int_{f_{j}}^{f_{j+1}}\\!{\rm d}N$ and the binary
parameters are generated randomly according to the distribution

$p(\mathcal{M}_{z},\eta,z,f)\propto\frac{{\rm d}N}{{\rm d}\mathcal{M}_{z}{\rm
d}\eta{\rm d}z{\rm d}\ln f}\,.$ (15)

This is the expression that was used to calculate the distributions of binary
parameters shown in Fig. 3 for a sampling with $p_{\rm BH}=0.17$ of $4\times
10^{5}$ SMBH binaries with redshifts $0<z<3$, mass ratios $\eta>0.01$, and GW
frequencies $1\,{\rm nHz}<f<30{\rm nHz}$. We display only binaries with chirp
mass $\mathcal{M}_{0}>10^{9}(f/{\rm nHz})^{-\frac{4}{3}}$, so that
$\bar{N}(f_{j})>\mathcal{O}(10^{4})$ for each $1\,{\rm nHz}$ frequency bin.

In order to generate the statistical distributions of $\Omega_{\rm 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(\Omega^{(1)},f)$ for individual sources $\Omega_{\rm 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\sigma$ and
$2\sigma$ 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\gtrsim 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\sigma$ CL range is due to the fact that the spectrum has
a long tail at high values of $\Omega_{\rm GW}$ generated by occasional nearby
binaries: the median spectrum always lies within the $1\sigma$ 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 $\Omega_{\rm GW}\propto f^{\frac{2}{3}}$ expectation in higher
frequency bins, while a few bins display sharp peaks.

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,
$\Omega^{(1)}_{\rm GW}\propto D_{L}^{-2}$, by Eqs. (10) and (8), and the
probability of finding an event at $D_{L}$ is $P^{(1)}(D_{L})\propto
D_{L}^{2}$. Thus $\langle\Omega^{(1)}_{\rm GW}\rangle\propto\int^{D_{L,\rm
max}}_{D_{L,\rm min}}{\rm d}D_{L}P(D_{L})\Omega^{(1)}_{\rm GW}\propto D_{L,\rm
max}$, where, in this simplified approach, $D_{L,\rm max}$ is some large
luminosity distance at which $\Omega^{(1)}_{\rm GW}$ gets suppressed, and
$D_{L,\rm min}$ is the distance to the nearest possible source. On the other
hand, $\langle(\Omega^{(1)}_{\rm 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)}(D_{L})$ that $\Omega^{(1)}_{\rm GW}$ has a
relatively flat power-law tail at large values,

$P^{(1)}(\Omega)\propto\Omega^{-\frac{5}{2}},\quad\mbox{when}\quad\Omega\to\infty\,.$
(16)

Since the closest distance to massive BH binaries is constrained, this tail
will be cut off at $\Omega^{(1)}_{\rm 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 $\Omega_{\rm
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\propto\Omega^{-\frac{5}{2}}$ (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. 8 in
Appendix A. In summary, even if a bin has a high number of contributing events
(i.e., $\bar{N}\to\infty$), the distribution of $\Omega_{\rm GW}$ does not
converge to a Gaussian because it retains its $\Omega_{\rm GW}^{-5/2}$ tail.
Since the variance is not well behaved, we define the width of the
distribution $\Delta_{\Omega_{\rm GW}}$ as the width of the 68% confidence
interval, as in Fig. 5. The width-to-mean ratio scales as:

$\Delta_{\Omega_{\rm GW}}/\left\langle\Omega_{\rm
GW}\right\rangle\propto\bar{N}^{-\frac{1}{3}}\,,$ (17)

when $\bar{N}\gg 1$, as is shown in Appendix A. This scaling is slower than
the typical $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
$\theta$ are

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

where

$|h(f)|=\frac{4\mathcal{M}_{z}^{\frac{5}{3}}}{D_{L}}(\pi f)^{\frac{2}{3}}$
(19)

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

$\displaystyle I(f)=|h_{L}(f)|^{2}+|h_{R}(f)|^{2}\,,$ (20) $\displaystyle
V(f)=|h_{L}(f)|^{2}-|h_{R}(f)|^{2}\,.$

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

$\mathcal{P}_{\rm circ}(f)=\frac{\sum_{i}V_{i}(f)}{\sum_{i}I_{i}(f)}\,,$ (21)

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

If the GW signal is dominated by a single source, then $\mathcal{P}_{\rm
circ}$ depends only on the inclination angle and large circular polarizations
are preferred, with $\mathcal{P}_{\rm circ}>0.87(0.2)$ at the 68% CL (95% CL).
On the other hand, if the GW signal is dominated by several ($N_{\rm
dom}\gtrsim 10$) sources of comparable strengths, then the $\mathcal{P}_{\rm
circ}(f)$ will be approximately Gaussian with a width determined by
fluctuations in $\sigma_{V}$. Since $\sigma_{V}/\langle
I\rangle_{\theta}=1.17$ for a single source, we can estimate that (see
Appendix A for details)

$\sigma_{\mathcal{P}_{\rm circ}}\approx\sigma_{V}/\langle
I\rangle_{\theta}\approx 1.17/\sqrt{N_{\rm dom}}\,.$ (22)

The lower-left panel of Fig. 5 illustrates the distributions of the circular
polarization $\mathcal{P}_{\rm circ}(f)$ found in the sample of 1600 Monte
Carlo realizations of the SMBH binary population with $M>10^{9}M_{\odot}$ and
$\eta>0.01$ whose frequency spectra were illustrated in Fig. 5. We see that
the mean value of $\mathcal{P}_{\rm circ}(f)\approx 0$, as expected, but large
statistical fluctuations are possible even at the $1\sigma$ 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
$\mathcal{P}_{\rm 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 $\sigma=0.13$ (blue) and $\sigma=0.23$ (orange). These fits are
very accurate as also demonstrated in the right panel of Fig. 8 in Appendix A.

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_{\rm BH}$ is not universal. By suppressing $p_{\rm BH}$, for instance,
at $M>10^{9}M_{\odot}$ and enhancing $p_{\rm BH}$ for
$10^{9}M_{\odot}>M>10^{6}M_{\odot}$, we can accommodate the IPTA common-
spectrum effect. In this case, the GW energy spectrum $\Omega_{\rm GW}$
becomes smoother and approaches the naive $\Omega_{\rm GW}\propto
f^{\frac{2}{3}}$ behavior while suppressing the circular polarization
$\mathcal{P}_{\rm circ}$. Moreover, in this case, $p_{\rm BH}$ for the lighter
binaries would be enhanced, which would increase the number of signals in the
sensitivity range of LISA and possibly AEDGE.

Figure 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 $\mathcal{M}_{0}$ and redshift $z$
(left panel) and of the chirp mass and symmetric mass ratio $\eta$ (right
panel). The merger rates assume the abundance of IMBHs described by Eq. (6)
for BHs heavier than $10^{3}M_{\odot}$ and a universal merger efficiency
factor of $p_{\rm BH}=0.17$ inferred from IPTA data, as discussed in the text.

###  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_{\rm insp}=\int_{\mathcal{T}+\tau_{\rm min}}^{\infty}\\!\\!{\rm
d}\tau\int{\rm d}\lambda\,p_{\rm det}\\!\left[\frac{{\rm SNR}_{c}}{{\rm
SNR}(\tau,z)}\right]\,,$ (23)

where we account only for binaries whose coalescence time $\tau$ is longer
than $\tau_{\rm min}=1\,$day and we estimate the SNR for a detector
characterized by the noise power spectrum $S_{n}(f)$ as

${\rm SNR}(\tau,z)=\sqrt{\int_{0}^{\mathcal{T}}{\rm
d}t\,\frac{2|h(\tau-t,z)|^{2}}{S_{n}(f(\tau-t,z))}}\,.$ (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_{\rm merg}=\mathcal{T}\int{\rm d}\lambda\,p_{\rm det}\\!\left[\frac{{\rm
SNR}_{c}}{{\rm SNR}(z)}\right]\,,$ (25)

where

${\rm SNR}(z)=\sqrt{\int_{f(\tau_{\rm min})}^{\infty}{\rm
d}f\,\frac{4|\tilde{h}(f)|^{2}}{S_{n}(f)}}$ (26)

The cut at $f(\tau_{\rm min})$ here implies that we only consider the SNR only
from the last day of the signal. In both cases, the detection probability
$p_{\rm 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 ${\rm SNR}_{c}=8$ for the detection
threshold and $\mathcal{T}=1\,{\rm year}$ for the observation time.

In the monochromatic limit, Eq. (24) is simplified to ${\rm
SNR}=|h(f)|\sqrt{2\mathcal{T}/S_{n}(f)}$. 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 ${\rm SNR}=8$. We see that IPTA30
can potentially resolve several, $\mathcal{O}(10)$, monochromatic GW signals
from SMBH binaries.

Figure 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 $\mathcal{M}_{0}$ and redshift
$z$ (upper panels) and of chirp mass and symmetric mass ratio $\eta$ (lower
panels). The merger rates assume the abundance of IMBHs described by Eq. (6)
for BHs heavier than $10^{3}M_{\odot}$ and a universal merger efficiency
factor of $p_{\rm BH}=0.17$ inferred from the IPTA data.

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
$10^{3}M_{\odot}$ is less than a day when they enter the AEDGE sensitivity
window.666In 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. 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
$\mathcal{M}_{0}=10^{4}M_{\odot}$ reflecting the binary population, that
increases towards lower masses, and the cut-off at $10^{3}M_{\odot}$.
Similarly, because the GW signal (10) from inspirals does not explicitly
depend on $\eta$, the $\eta$ 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\sim
4$.

We find that, at $z<10$, there are a total of $\mathcal{O}(10^{7})$ near-
monochromatic sources in the frequency range of $10^{-5}{\rm Hz}<f<0.1{\rm
Hz}$. 777We emphasize, however, that there could be lower-mass BHs that are
remnants of Population III stars, which are not included in our analysis. 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_{\rm BH}$ introduces additional
uncertainties in the merger rate at higher redshifts ($z\gtrsim 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,
$\mathcal{M}_{0}>10^{4}M_{\odot}$, while AEDGE can see lighter binaries,
$\mathcal{M}_{0}<2\times 10^{3}M_{\odot}$. LISA could (in principle) detect
even heavier mergers, namely, of $\mathcal{M}_{0}\lesssim 10^{8}M_{\odot}$,
but since such mergers are so rare it is unlikely that LISA will see any
mergers above $10^{6}M_{\odot}$. 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 (cf. Fig. 1) while
for LISA the $z$ distribution peaks at slightly lower $z$ because the SNR of
$\mathcal{M}_{0}\sim 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 AION-km and ET we find that the total expected
number of events per year is less than three.

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. 2022). 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 & Loeb
2022; Fragione et al. 2019; Fragione & Kocsis 2018) 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, $p_{\rm 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 mass-independent value of $p_{\rm
BH}$ to the PTA signal, finding $p_{\rm 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_{\odot}$ and about 10% from
masses $<10^{9}M_{\odot}$. The PTA mergers would have redshifts ${\cal O}(1)$
and have mass asymmetries $\gtrsim 10$.

The number of mergers contributing most of the PTA signal at frequencies
${\cal O}(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 $p_{\rm BH}$ as for the PTA
signal, there would be observable signals from mergers with total masses
$\in(10^{3},10^{6})M_{\odot}$ in the LISA experiment and from mergers with
total masses $\in(10^{3},10^{5})M_{\odot}$ in the AEDGE experiment. Data from
these experiments will be able to check the accuracy of the constancy of
$p_{\rm BH}$. In principle, their signals could be even larger than our
estimates if $p_{\rm 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_{\rm BH}$ and, hence, a lower event rate cannot
be excluded.

###### Acknowledgements.

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

  * Abbott et al. (2016) Abbott, B. P. et al. 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2021) Abbott, R. et al. 2021 [arXiv:2111.03634]
  * Aghanim et al. (2020) Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6, [Erratum: Astron.Astrophys. 652, C4 (2021)]
  * Ajith et al. (2008) Ajith, P. et al. 2008, Phys. Rev. D, 77, 104017, [Erratum: Phys.Rev.D 79, 129901 (2009)]
  * Akiyama et al. (2019) Akiyama, K. et al. 2019, Astrophys. J. Lett., 875, L1
  * Akiyama et al. (2022) Akiyama, K. et al. 2022, Astrophys. J. Lett., 930, L12
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P. et al. 2017 [arXiv:1702.00786]
  * Antoniadis et al. (2022) Antoniadis, J. et al. 2022, Mon. Not. Roy. Astron. Soc., 510, 4873
  * Arzoumanian et al. (2020) Arzoumanian, Z. et al. 2020, Astrophys. J. Lett., 905, L34
  * Arzoumanian et al. (2023) Arzoumanian, Z. et al. 2023 [arXiv:2301.03608]
  * Ashoorioon et al. (2022) Ashoorioon, A., Rezazadeh, K., & Rostami, A. 2022, Phys. Lett. B, 835, 137542
  * Badurina et al. (2021) Badurina, L., Buchmueller, O., Ellis, J., et al. 2021, Phil. Trans. A. Math. Phys. Eng. Sci., 380, 20210060
  * Badurina et al. (2020) Badurina, L. et al. 2020, JCAP, 05, 011
  * Barkana & Loeb (2001) Barkana, R. & Loeb, A. 2001, Phys. Rept., 349, 125
  * Bécsy et al. (2022) Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022, Astrophys. J., 941, 119
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307
  * Belgacem & Kamionkowski (2020) Belgacem, E. & Kamionkowski, M. 2020, Phys. Rev. D, 102, 023004
  * Benson et al. (2005) Benson, A. J., Kamionkowski, M., & Hassani, S. H. 2005, Mon. Not. Roy. Astron. Soc., 357, 847
  * Berti et al. (2006) Berti, E., Cardoso, V., & Will, C. M. 2006, Phys. Rev. D, 73, 064030
  * Blanco-Pillado et al. (2021) Blanco-Pillado, J. J., Olum, K. D., & Wachter, J. M. 2021, Phys. Rev. D, 103, 103512
  * Blas & Jenkins (2022) Blas, D. & Jenkins, A. C. 2022, Phys. Rev. Lett., 128, 101103
  * Blasi et al. (2021) Blasi, S., Brdar, V., & Schmitz, K. 2021, Phys. Rev. Lett., 126, 041305
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, Astrophys. J., 379, 440
  * Boroson & Lauer (2009) Boroson, T. A. & Lauer, T. R. 2009, Nature, 458, 53
  * Buchmuller et al. (2020) Buchmuller, W., Domcke, V., & Schmitz, K. 2020, Phys. Lett. B, 811, 135914
  * Chadayammuri et al. (2022) Chadayammuri, U., Bogdan, A., Ricarte, A., & Natarajan, P. 2022 [arXiv:2212.04693]
  * Chen et al. (2017) Chen, S., Middleton, H., Sesana, A., Del Pozzo, W., & Vecchio, A. 2017, Mon. Not. Roy. Astron. Soc., 468, 404, [Erratum: Mon.Not.Roy.Astron.Soc. 469, 2455–2456 (2017)]
  * Chen et al. (2021) Chen, S. et al. 2021, Mon. Not. Roy. Astron. Soc., 508, 4970
  * Conneely et al. (2019) Conneely, C., Jaffe, A. H., & Mingarelli, C. M. F. 2019, Mon. Not. Roy. Astron. Soc., 487, 562
  * De Luca et al. (2021) De Luca, V., Franciolini, G., & Riotto, A. 2021, Phys. Rev. Lett., 126, 041303
  * De Paolis et al. (2002) De Paolis, F., Ingrosso, G., & Nucita, A. A. 2002, Astron. Astrophys., 388, 470
  * De Paolis et al. (2004) De Paolis, F., Ingrosso, G., & Nucita, A. A. 2004, Astron. Astrophys., 426, 379
  * Dodelson (2003) Dodelson, S. 2003, Modern Cosmology (Amsterdam: Academic Press)
  * Eisenstein & Hu (1998) Eisenstein, D. J. & Hu, W. 1998, Astrophys. J., 496, 605
  * El-Neaj et al. (2020) El-Neaj, Y. A. et al. 2020, EPJ Quant. Technol., 7, 6
  * Ellis & Lewicki (2021) Ellis, J. & Lewicki, M. 2021, Phys. Rev. Lett., 126, 041304
  * Enoki et al. (2004) Enoki, M., Inoue, K. T., Nagashima, M., & Sugiyama, N. 2004, Astrophys. J., 615, 19
  * Enoki & Nagashima (2007) Enoki, M. & Nagashima, M. 2007, Prog. Theor. Phys., 117, 241
  * Erickcek et al. (2006) Erickcek, A. L., Kamionkowski, M., & Benson, A. J. 2006, Mon. Not. Roy. Astron. Soc., 371, 1992
  * Finn & Chernoff (1993) Finn, L. S. & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198
  * Fragione et al. (2019) Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, Mon. Not. Roy. Astron. Soc., 488, 47
  * Fragione & Kocsis (2018) Fragione, G. & Kocsis, B. 2018, Phys. Rev. Lett., 121, 161103
  * Fragione & Loeb (2022) Fragione, G. & Loeb, A. 2022 [arXiv:2212.04056]
  * Gerosa et al. (2019) Gerosa, D., Ma, S., Wong, K. W. K., et al. 2019, Phys. Rev. D, 99, 103004
  * Goncharov et al. (2021) Goncharov, B. et al. 2021, Astrophys. J. Lett., 917, L19
  * Graham et al. (2015) Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74
  * Graham et al. (2017) Graham, P. W., Hogan, J. M., Kasevich, M. A., Rajendran, S., & Romani, R. W. 2017 [arXiv:1711.02225]
  * Greene et al. (2020) Greene, J. E., Strader, J., & Ho, L. C. 2020, Ann. Rev. Astron. Astrophys., 58, 257
  * Hellings & Downs (1983) Hellings, R. W. & Downs, G. S. 1983, Astrophys. J. Lett., 265, L39
  * Hogan et al. (2011) Hogan, J. M. et al. 2011, Gen. Rel. Grav., 43, 1953
  * Holgado et al. (2018) Holgado, A. M., Sesana, A., Sandrinelli, A., et al. 2018, Mon. Not. Roy. Astron. Soc., 481, L74
  * Hotinli et al. (2019) Hotinli, S. C., Kamionkowski, M., & Jaffe, A. H. 2019, Open J. Astrophys., 2, 8
  * Hütsi et al. (2023) Hütsi, G., Raidal, M., Urrutia, J., Vaskonen, V., & Veermäe, H. 2023, Phys. Rev. D, 107, 043502
  * Iguchi et al. (2010) Iguchi, S., Okuda, T., & Sudou, H. 2010, Astrophys. J. Lett., 724, L166
  * Izquierdo-Villalba et al. (2021) Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2021, Mon. Not. Roy. Astron. Soc., 509, 3488
  * Janssen et al. (2015) Janssen, G. et al. 2015, PoS, AASKA14, 037
  * Kaiser & McWilliams (2021) Kaiser, A. R. & McWilliams, S. T. 2021, Class. Quant. Grav., 38, 055009
  * Kato & Soda (2016) Kato, R. & Soda, J. 2016, Phys. Rev. D, 93, 062003
  * Kelley et al. (2017a) Kelley, L. Z., Blecha, L., & Hernquist, L. 2017a, Mon. Not. Roy. Astron. Soc., 464, 3131
  * Kelley et al. (2017b) Kelley, L. Z., Blecha, L., Hernquist, L., Sesana, A., & Taylor, S. R. 2017b, Mon. Not. Roy. Astron. Soc., 471, 4508
  * King et al. (2013) King, O. G., Hovatta, T., Max-Moerbeck, W., et al. 2013, Mon. Not. Roy. Astron. Soc., 436, 114
  * Kocsis & Sesana (2011) Kocsis, B. & Sesana, A. 2011, Mon. Not. Roy. Astron. Soc., 411, 1467
  * Kocsis & Sesana (2011) Kocsis, B. & Sesana, A. 2011, MNRAS, 411, 1467
  * Kohri & Terada (2021) Kohri, K. & Terada, T. 2021, Phys. Lett. B, 813, 136040
  * Kovačević et al. (2022) Kovačević, A. B., Songsheng, Y.-Y., Wang, J.-M., & Popovic, L. C. 2022, Astron. Astrophys., 663, A99
  * Lacey & Cole (1993) Lacey, C. G. & Cole, S. 1993, Mon. Not. Roy. Astron. Soc., 262, 627
  * Lewicki & Vaskonen (2021) Lewicki, M. & Vaskonen, V. 2021, Eur. Phys. J. C, 83, 168
  * Lippai et al. (2009) Lippai, Z., Frei, Z., & Haiman, Z. 2009, ApJ, 701, 360
  * Merritt (2013) Merritt, D. 2013, Class. Quant. Grav., 30, 244005
  * Miller et al. (2015) Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, Astrophys. J., 799, 98
  * O’Neill et al. (2022) O’Neill, S. et al. 2022, Astrophys. J. Lett., 926, L35
  * Peters (1964) Peters, P. C. 1964, Phys. Rev., 136, B1224
  * Peters & Mathews (1963) Peters, P. C. & Mathews, J. 1963, Phys. Rev., 131, 435
  * Phinney (2001) Phinney, E. S. 2001 [arXiv:astro-ph/0108028]
  * Press & Schechter (1974) Press, W. H. & Schechter, P. 1974, Astrophys. J., 187, 425
  * Rieger & Mannheim (2000) Rieger, F. M. & Mannheim, K. 2000, Astron. Astrophys., 359, 948
  * Sathyaprakash et al. (2012) Sathyaprakash, B. et al. 2012, Class. Quant. Grav., 29, 124013, [Erratum: Class.Quant.Grav. 30, 079501 (2013)]
  * Sato-Polito & Kamionkowski (2022) Sato-Polito, G. & Kamionkowski, M. 2022, Phys. Rev. D, 106, 023004
  * Sesana et al. (2004) Sesana, A., Haardt, F., Madau, P., & Volonteri, M. 2004, Astrophys. J., 611, 623
  * Sesana et al. (2008) Sesana, A., Vecchio, A., & Colacino, C. N. 2008, Mon. Not. Roy. Astron. Soc., 390, 192
  * Sesana et al. (2009) Sesana, A., Vecchio, A., & Volonteri, M. 2009, Mon. Not. Roy. Astron. Soc., 394, 2255
  * Taylor et al. (2016) Taylor, S. R., Huerta, E. A., Gair, J. R., & McWilliams, S. T. 2016, Astrophys. J., 817, 70
  * Valbusa Dall’Armi et al. (2023) Valbusa Dall’Armi, L., Nishizawa, A., Ricciardone, A., & Matarrese, S. 2023 [arXiv:2301.08205]
  * Valtonen et al. (2008) Valtonen, M. J. et al. 2008, Nature, 452, 851
  * Vaskonen & Veermäe (2021) Vaskonen, V. & Veermäe, H. 2021, Phys. Rev. Lett., 126, 051303
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177
  * Volonteri et al. (2021) Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nature Rev. Phys., 3, 732
  * Wang et al. (2021) Wang, Y., Pardo, K., Chang, T.-C., & Doré, O. 2021, Phys. Rev. D, 103, 084007
  * Wyithe & Loeb (2003) Wyithe, J. S. B. & Loeb, A. 2003, Astrophys. J., 590, 691

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

###  A.1 Distribution of $\Omega_{\rm 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 $\Omega_{\rm GW}$ coming from a single
inspiralling binary:

$\Omega_{\rm GW}^{(1)}(\@vec{\theta})=\left|\frac{{\rm d}\tau}{{\rm d}\ln
f}\right|^{-1}\frac{2\pi}{5}\frac{f^{3}|\tilde{h}(f)|^{2}}{\rho_{c}}\,,$ (27)

where $\@vec{\theta}\equiv\\{\mathcal{M}_{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

$\Omega_{\rm GW}(f)\approx\sum_{i}\Omega_{{\rm
GW}}^{(1)}(\@vec{\theta}_{i})\delta(\ln f/f_{i})\,.$ (28)

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

$\displaystyle P^{(1)}(\Omega|f)$ $\displaystyle=\frac{1}{n(f)}\int{\rm
d}\lambda\left|\frac{{\rm d}\tau}{{\rm d}\ln
f}\right|\delta\left(\Omega-\Omega^{(1)}_{\rm GW}\right)$ (29)
$\displaystyle=\frac{\mathcal{F}(\Omega
f^{-\frac{10}{3}})}{f\Omega^{\frac{3}{2}}n(f)}\,,$

where we define the following function:

$\displaystyle\mathcal{F}(x)\equiv\int{\rm d}z$ $\displaystyle{\rm
d}\eta\bigg{[}\frac{\partial\lambda(\mathcal{M}_{z},z)}{\partial\mathcal{M}_{z}\partial\eta\partial
z}$ (30) $\displaystyle\
\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}}}$

and the spectral source density:

$n(f)\equiv\int{\rm d}\mathcal{M}_{z}{\rm d}\eta{\rm d}z\frac{{\rm d}N}{{\rm
d}\mathcal{M}_{z}{\rm d}\eta{\rm d}z{\rm d}\ln f}\,.$ (31)

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 $(f_{j},f_{j+1})$ can be
obtained from

$\displaystyle\Omega_{{\rm GW}}(f_{j})$
$\displaystyle\equiv\frac{1}{\ln(f_{j+1}/f_{j})}\int^{f_{j+1}}_{f_{j}}\Omega_{{\rm
GW}}(\@vec{\theta};f){\rm d}\ln f$ (32)
$\displaystyle=\frac{1}{\ln(f_{j+1}/f_{j})}\sum^{N(f_{j})}_{i=1}\Omega^{(1)}_{{\rm
GW},i}\,,$

where the number of binaries $N(f_{j})$ is drawn from a Poisson distribution
with the expected value,

$\bar{N}(f_{j})=\int_{f_{j}}^{f_{j+1}}\\!{\rm d}\ln f\,n(f)\,,$ (33)

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

The moment generating function of $\Omega_{{\rm GW}}$ is given by:

$\displaystyle M_{\Omega_{\rm GW}}(s)$
$\displaystyle\equiv\langle\exp\left(s\Omega_{\rm GW}\right)\rangle$ (34)
$\displaystyle=\sum_{N\geq
0}p_{N}\prod^{N}_{i=1}\left\langle\exp\left(s\Omega^{(1)}_{\rm
GW}\right)\right\rangle$ $\displaystyle=\sum_{N\geq
0}\frac{{\bar{N}}^{N}}{N!}e^{-\bar{N}}\left(M_{\Omega^{(1)}_{\rm
GW}}(s)\right)^{\bar{N}}$
$\displaystyle=\exp\left[\bar{N}\left(M_{\Omega^{(1)}_{\rm
GW}}(k)-1\right)\right]\,.$

Thus, the cumulant generating function of $\Omega_{{\rm GW}}(f_{j})$, that is,

$\displaystyle K_{\Omega_{\rm GW}}(s)$ $\displaystyle\equiv\ln M_{\Omega_{\rm
GW}}(s)\,$ (35) $\displaystyle=\int_{f_{j}}^{f_{j+1}}{\rm
d}N\left[\exp(s\,\Omega^{(1)}_{\rm GW})-1\right]\,.$

is proportional to the generating function of $M_{\Omega^{(1)}_{\rm GW}}(s)$,
implying that the $n$th cumulant of $\Omega_{{\rm GW}}$ is

$\displaystyle\kappa_{n}[\Omega_{\rm GW}]=\bar{N}\langle(\Omega^{(1)}_{\rm
GW})^{n}\rangle\,=\int_{f_{j}}^{f_{j+1}}{\rm d}N(\Omega^{(1)}_{\rm
GW})^{n}\,.$ (36)

Therefore, the mean $\kappa_{1}(\Omega_{\rm GW})\equiv\langle\Omega_{\rm
GW}\rangle$ matches the expectation value given in Eq. (8), while the
variance,

$\kappa_{2}(\Omega_{\rm GW})\equiv\langle(\delta\Omega_{\rm
GW})^{2}\rangle=\bar{N}\langle(\Omega^{(1)}_{\rm GW})^{2}\rangle$ (37)

is divergent due to the long $P^{(1)}(\Omega)\propto\Omega^{-\frac{5}{2}}$
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 $\Omega_{{\rm GW}}$, i.e., the higher cumulants are not
diminished when summing several instances of
$P^{(1)}(\Omega)\propto\Omega^{-\frac{5}{2}}$.

The moment-generating function $M_{\Omega_{\rm 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
$\Omega_{\rm GW}(f_{j})$ can be obtained by an inverse Laplace transform of
the moment generating function given in Eq. (34). It can be expressed as

$\displaystyle P(\Omega_{\rm GW})$
$\displaystyle=\frac{1}{2\pi\bar{N}}\int^{\infty}_{-\infty}{\rm
d}s\,e^{is\Omega_{\rm GW}/\bar{N}}\times$ (38)
$\displaystyle\times\exp\left(\bar{N}\left\langle\exp\left(-is\Omega^{(1)}_{\rm
GW}/\bar{N}\right)-1\right\rangle\right)\,,$

where the average is taken over the single event distribution $P^{(1)}$.

#### The large $N$ limit. –

The central limit theorem (and Eq. (37)) dictates that when $N\to\infty$, then
the distribution (32) should approach a Gaussian with its relative width
scaling as $\sigma_{\Omega_{\rm GW}}/\langle\Omega_{\rm 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)}(\Omega)\stackrel{{\scriptstyle\Omega\to\infty}}{{\sim}}C\Omega^{-\frac{5}{2}},$
(39)

where $C$ is a constant. The second moment $\langle\Omega^{2}\rangle$
diverges. We assume that the small $\Omega$ behavior is such that the first
moment $\bar{\Omega}^{(1)}\equiv\langle\Omega^{(1)}\rangle$ is finite. The
$\bar{N}\to\infty$ asymptotic of the average in Eq. (38) 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 (39):

$\displaystyle\bar{N}\left\langle\exp\left(-is\frac{\Omega^{(1)}}{\bar{N}}\right)-1\right\rangle$
$\displaystyle=$ $\displaystyle-
is\bar{\Omega}^{(1)}+\bar{N}\left\langle\exp\left(-is\frac{\Omega^{(1)}}{\bar{N}}\right)-1+is\frac{\Omega^{(1)}}{\bar{N}}\right\rangle$
$\displaystyle\stackrel{{\scriptstyle N\to\infty}}{{\sim}}$ $\displaystyle-
is\bar{\Omega}^{(1)}+e^{i\pi\frac{3}{4}}\frac{\sqrt{\pi}4C}{3\sqrt{\bar{N}}}s^{\frac{3}{2}}\,.$
(40)

In summary, the $N\to\infty$ distribution asymptotes to

$\displaystyle P(\Omega)$
$\displaystyle\approx(C\bar{N})^{-\frac{2}{3}}\mathcal{P}\left((C\bar{N})^{-\frac{2}{3}}(\Omega-\bar{\Omega})\right)\,,$
(41)

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

$\displaystyle\mathcal{P}\left(x\right)\equiv{\rm
Re}\frac{e^{i\pi/6}}{\pi}\int^{\infty}_{0}{\rm
d}s\,\exp\left(ie^{i\pi/6}xs-\frac{4\sqrt{\pi}}{3}s^{\frac{3}{2}}\right)$ (42)

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

  1. 1.

Since $\mathcal{P}(x)\sim x^{-5/2}$ when $x\to\infty$, the power-law behavior
of the tail is preserved in the large $\bar{N}$ limit:

$P(\Omega)\stackrel{{\scriptstyle\Omega\gg\bar{\Omega}}}{{\approx}}\frac{C\bar{N}}{(\Omega-\bar{\Omega})^{\frac{5}{2}}}\,.$
(43)

  2. 2.

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

$P^{(AN)}(\Omega)\\!=\\!A^{-\frac{2}{3}}P^{(N)}\left(\\!A^{-\frac{2}{3}}(\Omega\\!-\\!\bar{\Omega}^{(AN)})\\!+\\!\bar{\Omega}^{(N)}\\!\right).$
(44)

  3. 3.

Although the variance diverges, we can define the width of the distribution
$\Delta_{\Omega}\equiv\Omega_{+}-\Omega_{-}$ from the confidence interval
$(\Omega_{-},\Omega_{+})$ centred around, e.g., the median. Eq. (44) then
implies the scaling

$\Delta_{\Omega}\propto\bar{N}^{2/3}\,.$ (45)

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

  4. 4.

Using the method of steepest decent, we find that

$\mathcal{P}\left(x\right)\stackrel{{\scriptstyle\Omega\ll\bar{\Omega}}}{{\approx}}\frac{\sqrt{\bar{\Omega}-\Omega}}{2\pi
C\bar{N}}e^{-\frac{(\bar{\Omega}-\Omega)^{3}}{12\pi(C\bar{N})^{2}}}\,,$ (46)

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)}(\Omega)\propto\Omega^{-\frac{5}{2}}\theta(\Omega-\Omega_{\min})\,,$
(47)

which mimics the statistics of $\Omega_{\rm GW}$ from a realistic BH
population found in section 3.2. In Fig. 8, we show the resulting
distributions of $\Omega/\bar{\Omega}$ if the signal consists of 10, 100, and
$10^{5}$ sources. The mean is given by $\bar{\Omega}\equiv 3N\Omega_{\rm
min}$. The distribution is computed by generating random realizations of the
source populations (shown by points) and by using (38) (shown by solid lines).
These two approaches are in excellent agreement. As expected, one can observe
that the total signal will inherit the long $\Omega^{-5/2}$ tail from the
distribution of individual sources (shown by the dashed line). At large
$\Omega$, the tail approaches Eq. (43) in all cases, while the asymptotic (41)
works well for $N=10^{5}$, but predicts a slightly too wide peak for $N\leq
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. 8 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.

Figure 8: Distributions of the total signal strengths from 1, 10, 100, and
$10^{5}$ sources whose independent signal strengths have identical
$\Omega^{-\frac{5}{2}}$ tails (left). Distribution of $\mathcal{P}_{\rm 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 $\mathcal{P}_{\rm circ}$

Considering now the distribution of the circular polarization

$\mathcal{P}_{\rm
circ}\equiv\frac{\sum^{N}_{i=1}V_{i}}{\sum^{N}_{i=1}I_{i}}\,,$ (48)

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.
(18-LABEL:eq:I,V). 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 $\langle\mathcal{P}_{\rm
circ}\rangle=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

$\mathcal{P}_{\rm
circ}=\frac{4c_{\theta}(1+c_{\theta}^{2})}{1+6c_{\theta}^{2}+c_{\theta}^{4}}\,,$
(49)

where $c_{\theta}\equiv\cos(\theta)$. The corresponding probability
distribution is shown as the dashed curve in Fig. 8 and can be seen to be
dominated by the region $|\mathcal{P}_{\rm circ}|\approx 1$.

If the signal is dominated by several sources, the analytic treatment of
$\mathcal{P}_{\rm circ}$ is complicated due to correlations between $V$ and
$I$. The mean and variance of $I$ and $V$ are given by

$\displaystyle\langle I\rangle$
$\displaystyle=\frac{4}{5}\langle|h|^{2}\rangle\,,\qquad\qquad\ \langle
V\rangle=0\,,$ (50) $\displaystyle\langle\delta I^{2}\rangle$
$\displaystyle=\frac{284}{315}\langle(\delta|h|^{2})^{2}\rangle\,,\ \
\langle\delta V^{2}\rangle=\frac{92}{105}\langle(\delta|h|^{2})^{2}\rangle\,,$

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 (37). The
latter is expected to be large. To estimate the fluctuations in
$\mathcal{P}_{\rm circ}$ in the limit when the number of sources is large, we
first assume that $I$ can be replaced by its average over inclinations
$\langle I\rangle_{\theta}\equiv(4/5)\sum_{i}|h_{i}|^{2}$, but not over signal
strengths, so that

$\mathcal{P}_{\rm circ}\approx V/\langle
I\rangle_{\theta}\equiv\sum^{N}_{i=1}p_{i}f_{V}(c_{\theta,i})\,,$ (51)

where $p_{i}\equiv(4/5)|h_{i}|^{2}/\langle I\rangle_{\theta}\in[0,1]$ is the
fractional contribution to the signal from source $i$ and
$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 $\langle I\rangle_{\theta}$ in the random fractions
$p_{i}$, which vary only in the range $[0,1]$. We find:

$\langle\delta\mathcal{P}_{\rm
circ}^{2}\rangle=\frac{115}{84}N\langle|p_{i}|^{2}\rangle\,.$ (52)

We note that $N\langle|p_{i}|^{2}\rangle$ takes values in the range $[1/N,1]$
and could be interpreted as a measure of the effective number of dominant
sources $1/N_{\rm 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

$\sigma_{\mathcal{P}_{\rm circ}}=1.17/\sqrt{N}\,,$ (53)

as also stated in Eq. (22). In the right panel of Fig. 8, approximate Gaussian
distributions with a width given by (53) and shown as solid lines are compared
to distributions obtained from explicitly generated random populations of
$\mathcal{P}_{\rm 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 (51) above can be improved by
expanding $\mathcal{P}_{\rm circ}$ in $\delta I\equiv I-\langle
I\rangle_{\theta}$, that is, $\mathcal{P}_{\rm circ}=V/\langle
I\rangle_{\theta}(1-\delta I/\langle I\rangle_{\theta}+\ldots)$. This allows
us to account for correlations between $V$ and $\delta I$, for instance, when
computing the variance (52),

$\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$
(54)

The corrections arising from higher powers of $\delta I/\langle
I\rangle_{\theta}$ are suppressed by increasing powers of $1/N$.