Probing primordial black holes at high redshift with future gravitational wave detectors

ABSTRACT

We analyse the detection prospects for potential Primordial Black Hole Binary (PBHB) populations buried in the Stellar-Origin Black Hole Binary (SOBHB) population inferred by the LVK collaboration. We consider different PBHB population scenarios and several future Gravitational Wave (GW) detectors. To separate the PBHB component from the SOBHB one, we exploit the prediction that the PBHB merger rate does not decline as fast as the SOBHB one at high redshift. However, only a tiny fraction of PBHB events may be resolved individually, and the sub-threshold events may yield an undetectable Stochastic GW Background (SGWB). For this reason, we determine the statistical significance of the PBHB contributions in the number of resolvable events seen in future Earth-based detectors and the SGWB measured at LISA. We quantify them in the limit that SOBHB population uncertainties are small, as one may optimistically expect at the time that future detectors will operate. In general, we find the synergy between these probes will consistently help assess whether or not a sizeable PBHB population is present.

1 INTRODUCTION

LIGO’s first Gravitational Wave (GW) detection labelled GW150914 (Abbott et al. 2016a, b) opened the gates to the world of GW astronomy. Since that detection, the sensitivity of LIGO detectors has increased considerably, and the Virgo (Acernese et al. 2023) and KAGRA (Akutsu et al. 2020) experiments have joined the network. At present, the LIGO-Virgo-KAGRA (LVK) network has identified over ninety GW events involving Black Holes (BHs) and Neutron Stars (NSs) (Abbott et al. 2023a) and has started constraining the statistical properties of the Stellar Mass Black Hole Binaries (SMBHB) population, although some bounds remain loose. Primarily, these limitations arise from the current detector sensitivity, which has allowed us to measure only a few events with high precision. The determination of the merger rate distribution at high redshift (z ≳ 1) is a striking example of these limitations. As a matter of fact, the properties of the population at z ≳ 1 are, currently, only guessed by using phenomenological models following the Star Formation Rate (SFR) (Madau & Dickinson 2014; Madau & Fragos 2017; Mangiagli et al. 2019), leaving open space for the presence of both long time delays in the BHB formation (Fishbach & Kalogera 2021) and a variety of populations with different redshift behaviours (Hütsi et al. 2021; Franciolini et al. 2022b; Antonelli et al. 2023) (see also refs. Cusin et al. 2020; Abbott et al. 2021a, 2023b; Périgois et al. 2021; Babak et al. 2023). While this might be a reasonable assumption if all the events observed by current detectors are BHBs of Stellar Origin (SOBHBs), at least in principle, different scenarios are possible, leading to different high redshift behaviours.

An intriguing alternative for BH formation is the possibility for BHs to form due to some cosmological processes occurring in the early Universe. These objects are typically dubbed Primordial BHs (PBHs) to differentiate from BHs produced in some late-time astrophysical processes. PBHs might form when strong scalar perturbations, e.g. generated by some inflationary model violating slow-roll, re-enter the Universe horizon, leading to the collapse of some regions of space (Zel’dovich & Novikov 1967; Hawking 1971; Carr & Hawking 1974; Carr 1975). Such a mechanism would make PBHs, and, in particular, PBH Binaries (PBHBs), totally independent of star formation processes. Different inflationary set-ups and early-universe histories can result, e.g. in vastly different PBH abundances and mass distributions (see e.g. refs. Carr et al. 2010; Khlopov 2010; García-Bellido 2017; Sasaki et al. 2018; Carr et al. 2021; Bagui et al. 2023 for reviews of PBH formation and constraints). Interestingly, depending on the formation scenario, PBHs can account for a substantial portion of the observed dark matter (DM) abundance (Carr & Kuhnel 2020; Carr et al. 2021; Villanueva-Domingo, Mena & Palomares-Ruiz 2021; Bagui et al. 2023) with experimental upper bounds depending on the formation mechanism and the PBH mass range. Furthermore, at least one of the two BHs in some SMBHBs might be of primordial origin, and, more in general, a PBHB population might contribute¹ to the events currently observed by LVK detector (Bird et al. 2016; Ali-Haïmoud et al. 2017; Clesse & García-Bellido 2017; Raidal et al. 2019; Hall, Gow & Byrnes 2020). In such a case, some statistical properties and observable signatures might radically differ from those arising when all SMBHBs are SOBHBs.

Despite their radically different origin and phenomenology, PBHs are elusive at current GW detectors. The challenge is partly rooted in the lack of unquestionable, discriminating predictions at low redshift [see ref. Franciolini et al. (2022a) for a systematic procedure to assess the origin of BHBs]. The predictions on PBHs and SOBHBs at low redshift are, indeed, loose due to the plethora of viable inflationary mechanisms and the numerous unknowns on the stellar-origin formation channels. On the contrary, a robust model-independent discrimination criterion exists at high redshift: at distances beyond the SFR peak (z ≃ 2), the SOBHB merger rate must fast decline, whereas the PBHB merger rate can keep growing (Ali-Haïmoud et al. 2017; Raidal et al. 2019; Young & Byrnes 2020; Hütsi et al. 2021; Atal et al. 2022; Bavera et al. 2022). Remarkably, while the SFR peak is beyond the reach of the present LVK interferometers, it will be in the range² of future Earth-based detectors (Koushiappas & Loeb 2017; Barsotti, McCuller & Fritschel 2018; Maggiore et al. 2020; De Luca et al. 2021; Ng et al. 2021; Martinelli et al. 2022; Branchesi et al. 2023; Franciolini et al. 2023; Ng et al. 2023). In addition, the Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017) might look for imprints of PBHs in the milli-Hertz band from individual events (Guo, Shu & Zhao 2019; Auclair et al. 2023; Bagui et al. 2023) and the Stochastic GW Background (SGWB)(Garcia-Bellido, Peloso & Unal 2017; Bartolo et al. 2019; Cai, Pi & Sasaki 2019; Unal 2019; Bagui et al. 2023). Moreover, LISA will be sensitive to the SGWB due to the incoherent superposition of the weak signals from the SMBHB population, which, including the contribution of binaries at high-redshift, brings information on the behaviour of the population above the SFR peak (Chen, Huang & Huang 2019; Cusin et al. 2020; Périgois et al. 2021; Bavera et al. 2022; Babak et al. 2023; Lewicki & Vaskonen 2023).

In this paper, we discuss the detection prospects for the PBHB population using its high-redshift behaviour (Martinelli et al. 2022; Ng et al. 2022a, b, 2023). Specifically, we study how future detectors might be able to identify PBHB populations beyond a certain Fiducial Population of SOBHB (based on ref. Babak et al. 2023), which is broadly compatible with the current LVK observations (Abbott et al. 2023a, 2023b). For this purpose, we consider some well-established PBHB population models (Atal et al. 2022; Bavera et al. 2022) and show how different detectors will complementarily probe their parameter spaces. We expect our qualitative results to be independent of our fiducial population and PBHB models. However, the quantitative outcomes are contingent upon our population selections, so it will be worth repeating our analysis when the statistical properties of the SMBHB populations become less uncertain than they are today.

The paper is structured as follows: in Section 2, we describe the PBHB (sub)population models, and in Section 3, we describe the methodology that we adopt. Our results are presented in Section 4, where we discuss the detectability of different PBHB populations with future Earth- and space-borne GW detectors. In particular, we demonstrate that the SGWB detection will give important complementary information on the presence of deviations from the fiducial model. We devote Section 5 to our final remarks and conclusions. Some technical details relevant to our analysis are described in Appendixes  A,  B,  D, and  E.

2 POPULATION MODELS

In this section, we present the population models that we analyse in this work. Following refs. Abbott et al. (2023b) and Babak et al. (2023), for a given population, the number of expected sources in a given interval of redshift and parameter space is given by

where |$\rm T^{Det}_{\rm Obs}$| is the detector observation time, dV_c/dz(z) is the differential comoving volume, R(z) is the merger rate, and p(θ|ξ) is the probability distribution function (PDF) for the source to have some specific values for the binary parameters (collectively denoted with θ) given some population hyperparameters (collectively denoted with ξ).³ Notice that since p(θ|ξ) is normalized, the number of events in a redshift interval [z_m, z_M] is given by

More in detail, the PDF term p(θ|ξ) can be expressed as

where m₁ and m₂ are the two masses, p(m₁, m₂|ξ_Mass) is the mass function depending on some hyperparameters ξ_Mass, p(θ_Angles|ξ_Angles) is the angle PDF, and p(θ_Spins|ξ_Spins) is the spin PDF. The sources are assumed to be isotropic in the sky, with inclination and polarization uniformly distributed in their considered prior. Eccentricity in the orbit is neglected throughout this work. The spin PDF is further expanded as

where |$p(a_1, a_2|\xi _{\rm Spin \, Amplitude})$| and |$p(\cos (t_1),\cos (t_2)|\xi _{\rm Spin \, Tilt})$| are the spin amplitude and spin tilt PDFs, with a_i and t_i (with i = 1, 2) denoting the normalized spin amplitude and the angle between the binary angular momentum and the spin of the body i, respectively.

In our analysis, the SMBHB population consists of the sum of the SOBHB and PBHB populations. We model each population using the framework outlined in equation (1) and neglect binaries with mixed origins. For the SOBHB population, we consider a fiducial population model, utilizing PDFs derived from the most recent LVK studies. We set their hyperparameters to the best-fitting values as determined by these studies (Abbott et al. 2023b). We remind that the LVK data, which we refer to as GWTC-3, provide no direct constraint at z ≳ 1. To extend our analysis to the redshift range 1 ≲ z ≲ 10, we assume that the SOBHB merger rate R_SOBHB tracks the SFR and that the PDFs remain redshift-independent. Specifically, we employ a phenomenological merger rate that closely follows the Madau-Dickinson SFR (Madau & Dickinson 2014) with negligible time delay (see Dominik et al. 2012; Mapelli et al. 2017; Neijssel et al. 2019). Appendix  A provides further details about our approach.

Due to their early-universe origin, PBHBs have different statistical properties than SOBHBs, and, in particular, the PBHB merger rate is not expected to track the SFR. While several possibilities exist in the literature (see e.g. ref. Atal et al. 2022), for the PBHB merger rate we adopt (Raidal et al. 2019; Young & Byrnes 2020)

which corresponds to a power law in cosmic time t(z) (Hütsi et al. 2021; Bavera et al. 2022), normalized so that in z = 0 it is ε-times smaller than the fiducial SOBHB merger rate R₀ at the same redshift. Alternatively, one can parametrize the rate by relating it to f_PBH ≡ Ω_PBH/Ω_DM, that is, the total PBH energy density over the DM energy density; see Appendix  B for details.

Fig. 1 shows the fiducial SOBHB merger rate and the PBHB merger rate for some values of ε. As long as ε ≪ 1, the PBHB population is subdominant within the LVK horizon (z ≲ 1). Thus, (up to a few outliers) the whole SMBHB population of GWTC-3 exhibits the SOBHB properties. On the other hand, the PBHB population may become relevant after the peak of the SFR. This is particularly evident from the dashed lines in the right panel of Fig. 1, which shows |$\hat{N}^{\rm Pop}_{z}$| number of events redshift z predicted by population Pop, which is given by

where |$\Delta N^{\rm Pop}_{z_{m}, z_{M}}$| is defined in equation (2) with the additional superscript specifying the considered population. By normalizing with Δz = z_M − z_m, we ensure that, for sufficiently small bin sizes, |$\hat{N}^{\rm Pop}_{z}$| is independent of the specific binning scheme used in the analysis.

For the PBHB mass distribution, we assume each of the two masses in the binary to be drawn from the same PDF,⁴ i.e. P(m₁) = P(m₂), which we set to be a Log-Normal (LN)

where μ_LN and σ_LN are the hyperparameters setting the position and width of the peak.

Finally, the PDFs of the angular variables and spins depend on the specific PBH formation mechanism and pairing (De Luca et al. 2019; Flores & Kusenko 2021; Saito et al. 2023). For concreteness, we consider a reasonable and practical option: both PBHB spin and angular distributions resemble those of the fiducial SOBHB population. This is also a pessimistic scenario, as reconstructing such PDFs would not help disentangle the two populations. However, the results of our analysis mildly depend on this choice, as discussed in Section 4.3.

3 METHODOLOGY

This section describes our methodology to identify a PBHB population component on top of the fiducial population consistent with GWTC-3. Our analysis studies the detectability of individual sources and SGWB with future GW detectors. Specifically, we consider LIGO |$\rm A^+$| and ET, with either 1 or 10 yr of observations, for the measurement of individual sources, and LISA, assuming either 4 or 10 yr of observations, for the SGWB detection and characterization. While our analysis concerns the PBHB population models presented in Section 2, a similar methodology could be applied to other PBH scenarios and GW detectors. For details on the detector characteristics and our numerical codes, see Section  D and the repository (Marcoccia 2023a), respectively.

3.1 Resolvable sources analysis

Our analysis of individually resolvable sources relies on measurements performed with Earth-based detectors. In particular, we check whether the presence of the PBHB population increases the expected number of detectable sources by more than |$3\, \sigma$| beyond the number predicted by the fiducial population. For this purpose, hereafter we define an event as ‘detectable’ when its signal-to-noise ratio (SNR) (see definition in equation 10) is larger than some threshold value SNR_Thr, which we set to 8.

Let us start by briefly reviewing the methods to evaluate the SNR associated with a GW signal. The signal, h, measured by any GW detector, is expressed as

where |$h^{+}_{ij}$| and |$h^{\times }_{ij}$| denote the two GW polarization modes, while |$F^{+}_{ij}$| and |$F^{\times }_{ij}$| are the detector pattern functions (for details, see e.g. ref. Maggiore 2000). The two GW modes can be further expanded as a combination of polarization tensors |$\rm {e}^{+}_{ij}$|⁠, |$\rm {e}^{\times }_{ij}$|⁠, and a waveform depending on the source parameters. In our study, we employ the IMRPhenomXHM waveform (Pratten et al. 2021), a phenomenological waveform from the IMRPhenom family (Ajith et al. 2007, 2008; Santamaria et al. 2010; Husa et al. 2016; Khan et al. 2016), offering a good compromise between quality and computation speed. With this choice, the signal depends on the parameters

where m₁ and m₂ are the masses of the two BHs in the detector frame, d_L is the luminosity distance, ϕ₀ is the binary’s initial phase, τ_c is the coalescence time, θ and ϕ are the binary’s latitude and longitude, ι is the angle between the binary’s angular momentum and the line of sight, ψ is the orientation, and χ₁ and χ₂ are the two (dimensionless) spin amplitudes projected on the orbital plane. Finally, the SNR for a given source and a given detector is defined as (Cutler & Flanagan 1994)

where S_n(f) is the detector strain sensitivity (see Section  D), while f_m and f_M are the minimal and maximal detector frequencies. Notice that in this equation, h depends on all the parameters listed in equation (9), but for a large sample of sources, only their average values matter in our analyses, at least at the leading order. This is why in several evaluations, e.g. the analytical estimates, we can average over both the spins (χ₁, χ₂) and angular (θ, ϕ, ι, ψ) variables; we dub this approximated SNR as |$\rm SNR_{\rm avg}$|⁠. However, as a check of robustness, in a few cases, we test our averaged-based results with those obtained without the average approximation, and we indicate the result of this precise evaluation simply as SNR (i.e. without any subscript).

As a first step, to fast probe the PBHB detectability in a broad part of their parameter space, we perform a semi-analytical analysis. For this purpose, we modify equation (6) by including a selection effect. Specifically, we define the expected number of resolvable sources at redshift z for the population Pop as

where, for any given detector, the selection function |$\theta _{\rm SNR_{\rm avg}}(\xi)$| is a Heaviside Θ function filtering the sources with |$\rm SNR_{\rm avg}$| larger than |$\rm SNR_{\rm Thr}=8$|⁠. The integrals in equation (11) are carried out numerically. Notice that due to the presence of the selection function, the integration over the ξ variables has to be computed explicitly.

We use equation (11) to set our (analytic) criterion for the identification of the PBHB component via resolvable sources at Earth-based detectors. We define the PBH contribution to be visible if, in a given bin in z, it satisfies the condition

In other words, the PBHB component can be separated from the fiducial SOBHB component if there exists a bin in z, in which the number of resolvable sources, |$\hat{N}^{\rm Res, PBHB}_{z }$|⁠, exceeds the number of detectable SOBHB sources, |$\hat{N}^{\rm Res, SOBHB}_{z}$|⁠, by 3|$\, \sigma$|⁠. In our case, the error comes from a Poissonian distribution, and this is why we have |$\Delta ^{\rm Res, Fid}_{z } = \sqrt{\hat{N}^{\rm Res, Fid}_{z}}$|⁠.

The semi-analytic analysis is fast but disregards two effects: the populations’ realization dependence and the impact of the angular and spin variables on the SNR. We quantify these effects by running a more sophisticated analysis on some PBHB population benchmarks. Specifically, we use the code in ref. Marcoccia (2023b) to sample over the PDF in equation (1) and generate n_r = 100 catalogues of the SOBHB fiducial population and one catalogue per PBHB benchmark population. Then, for every merger event predicted in the catalogues, we compute |$\rm SNR_{\rm avg}$| and the exact SNR. We define as |$\Delta \hat{\mathcal {N}}^{\rm Res, Pop}_{z_{\rm m}, z_{\rm M} }$| the number of events with |$\rm SNR_{\rm avg}\gt SNR_{Thr}$| in a given redshift interval z_m ≤ z ≤ z_M for the catalogue i of the population Pop. Similarly, we define as |$\Delta \overline{\mathcal {N}}^{\rm Res, Pop}_{z_{\rm m}, z_{\rm M, i} }$| the analogous quantity obtained with the detectability condition |$\rm SNR \gt SNR_{Thr}$|⁠. For convenience, we also introduce

The mismatches between |$\hat{\mathcal {N}}^{\rm Res, Pop}_{z, i}$| and |$\hat{N}_{z}^{\rm Res, Pop}$| highlight the effect of the realization dependence that our semi-analytic results neglect. However, we expect the mismatch to be statistically within the Poisson deviation from the mean, i.e.

where

Here, the index i runs over n, the number of realizations of each population scenario. Since we produce multiple realizations only of our fiducial population (namely n = n_r = 100), we focus on the realizations of this population to test equation (14). This allows us to prove |$\sigma ^{\rm Res, Pop}_z\simeq \Delta ^{\rm Res, Pop}_z$| and, in turn, to use |$\Delta ^{\rm Res, Pop}_z$| as a proxy of the realization dependences in our PBHB benchmarks.

The quantities |$\overline{\mathcal {N}}^{\rm Res, Pop}_{z_{\rm m}, z_{\rm M}, i }$| are useful to investigate the impact of the approximation |$\rm SNR_{\rm avg}$|⁠, in which the SNR is computed by averaging over the angles and spin. For this purpose, we calculate the quantities |$\overline{\mu }^{\rm \rm Res, Pop}_z$| and |$\overline{\sigma }^{\rm Pop, Res}_z$|⁠, given as in equations (15) and (16) but with the hat symbol replaced by the bar one. In the parameter regions where the approximation is satisfactory, |$\hat{\mu }^{\rm Res, Pop}_z$| and |$\overline{\mu }^{\rm Res, Pop}_z$| are expected to be practically equal.⁵ All these quantities can be computed for different values of |$\rm T^{\rm Det}_{\rm Obs}$| and several detector sensitivities. For concreteness, we consider LIGO A⁺ and ET, assuming |$T^{\rm Det}_{\rm Obs} =1, 10$|  yr of data. Moreover, to simplify the notation, hereafter we drop all the z subscripts in all these quantities and refer to the quantity, say q, as q^{Res, Pop}.

Let us comment on the assumptions of our semi-analytic analysis. The most relevant assumption pertains to the procedure to evaluate the right-hand side of equation (12). In particular, while a consistent analysis should account for both the realization error and the uncertainty in the model parameters, we evaluate equation (12) using the central values for all the fiducial population parameters without including their uncertainties. There are three main reasons behind this choice:

• The main message of this work is to stress the synergy between Earth-based and space-based GW detectors for what concerns assessing the presence of populations of high-redshift SMBHBs beyond our fiducial population. Thus, including further uncertainties in the analysis will quantitatively affect our results, but it will not change the message of this work.

• Current GW detections have only probed the Universe at z ≲ 1 so that the peak in the SOBHB population directly descends from imposing the population to follow the SFR at high redshift. Moreover, we have no information on the possible presence of time delays, which might shift the SOBHB peak position. All these uncertainties should also be included to perform a consistent analysis.

• With more measurements to come in the next few years (with improved sensitivity and possibly with more detectors joining the existing network), the determination of the model parameters will improve significantly. Currently, we have no reliable, precise estimate on the errors in the measured values of the fiducial population parameters at the end of the next LVK runs.⁶

For these reasons, we restrict ourselves to the case where the main uncertainty on |$\hat{N}^{\rm Res, Pop}$| comes from the Poissonian error. However, we show that this treatment is reasonably good by performing a more computationally demanding analysis on some benchmarks. We postpone the exhaustive treatment of the realization dependencies and population uncertainties to the time when more LVK data will be available.

As a final comment, we stress that the results presented in this work are sensitive to uncertainties on the astrophysical populations. For example, extra formation channels not captured in a time-delayed SFR might induce additional merger rate excesses and, in turn, some (likely small) systematics. However, depending on how different the mass/redshift PDFs are compared to the ones considered in this work, these channels could still be identified by more sophisticated analyses (De Luca et al. 2021; Bavera et al. 2022; Franciolini et al. 2022a, b; Ng et al. 2022b). For this reason, our approach constitutes a fast identification of PBH hints but does not substitute other methods that, at a later stage, can be tuned to confirm or disconfirm the found hints.

3.2 SGWB analysis

The SGWB from the fiducial SOBHB population, and its variation due to the PBHB contribution, is evaluated using the analytical approach summarized in Appendix  E and detailed in ref. Phinney (2001). It turns out that, within the LISA frequency band, the SGWB signal sourced by each population can be parametrized as a simple power law:

where β = 2/3 is the tilt of the GW power spectrum and α_Pop is the (log₁₀) of the amplitude at a reference (irrelevant) pivot frequency f_*. This result assumes each frequency bin to be highly populated by the GW signals due to binaries in circular orbits, with negligible environmental effects, and emitting GWs only (Cusin et al. 2020; Périgois et al. 2021; Bavera et al. 2022; Babak et al. 2024). Dropping any of these assumptions might induce modifications from the power-law behaviour.⁷ Since we are interested in evaluating the SGWB in the LISA frequency band |$f \in \lbrace 3 \times 10^{-5} \, , 0.5 \rbrace$| Hz, it is convenient to set f_* = 0.01 Hz.

Our strategy to probe the presence of the PBHB population via the SGWB measurement consists of reconstructing the overall signal with the above power-law template, and then checking whether the posteriors of the template parameters are compatible with those predicted by the fiducial SOBHB population up to some confidence level. As the aforementioned assumptions imply β = 2/3 for both SOBHB and PBH populations, only deviations from the SOBHB amplitude parameter are relevant to our goal. We calculate the reference SOBHB population amplitude by using equations (E5) and (E6) where, for practical purposes, the integral over z is cut off above a maximal redshift. Specifically, for the SOBHB population, we integrate up to z ≈ 10 and obtain α_SOBHB ≃ −12.02 at f_* = 0.01 Hz. With such a maximal-redshift choice, the estimate of the amplitude is accurate to within |$1 {{\ \rm per\ cent}}$| error (Babak et al. 2023). We follow the same procedure to compute α_PBH predicted by a PBH population with a given set of hyperparameter values. However, the (unresolved) PBHB signals do not die off as fast as the SOBHB ones. For this contribution, we set the cut-off at z = 100 corresponding to a |$\sim 10{{\ \rm per\ cent}}$| accuracy in the SGWB evaluation.⁸

To forecast the LISA posteriors on the power-law template parameters, we perform an analysis based on the Fisher Information Matrix (FIM) formalism. Given some data |$\tilde{d}(f)$|⁠, containing signal |$\tilde{s}(f)$| and noise |$\tilde{n}(f)$|⁠, which we assume to be Gaussian, with zero means, and characterized only by their variances,⁹ the (log-)likelihood can be written as

where f_min and f_max are the minimal and maximal frequencies measured by the detector, T is the total observation time, and |$D(f, \vec{\theta })$| is the model for the variance of the data, depending on some (signal and noise) parameters |$\vec{\theta }$|⁠. The best-fitting parameters |$\vec{\theta }_0$| are defined to maximize |$\log \mathcal {L}$|⁠:

which is clearly solved by |$D(f, \vec{\theta }_0) = \tilde{d}(f) \tilde{d}^{*}(f)$|⁠. Then, the FIM F_αβ is given by

By definition, the FIM F_αβ is the inverse of the covariance matrix C_αβ. As a consequence, estimates of the errors on the model parameters |$\vec{\theta }$| are obtained by computing |$\sqrt{ \rm {diag}(C_{\alpha \beta }) } = \sqrt{ \rm {diag}(F_{\alpha \beta }^{-1}) }$|⁠. Notice that LISA will measure three data streams. Under some simplifying assumptions, these data streams in the AET TDI basis are independent (see Section  D), and therefore the total Fisher matrix is given by

For what concerns the observation time T, we assume |$100 {{\ \rm per\ cent}}$| efficiency and impose |$T^{ \rm LISA}_{\rm Obs} = 4$| yr. For reference, we also show how results improve if the mission lifetime is extended to 10 yr. Let us assume that the data are expressed in Ω units, and we have factored the detector response out (for details, see Section  D). Then, |$D(f, \vec{\theta })$| can be expanded as

where |$\Omega _{\rm GW}(f, \vec{\theta }_s)$| is the template of the signal as a function of the frequency and the parameters |$\vec{\theta }_s=\lbrace \beta , \alpha _{\rm Tot}\rbrace$|⁠, and |$\Omega _{n}(f, \vec{\theta }_n)$| is the noise model as a function of the frequency and the noise parameters |$\vec{\theta }_n=\lbrace A,P \rbrace$|⁠. For what concerns the noise, we use the analytical two-parameter model commonly used in the literature (for details, see Section  D).

In our case, the signal is a sum of two contributions, one for the fiducial population and one for the PBHB population, each described by the template in equation (17). Given the complete degeneracy between these two components, we cannot measure them independently, but only the overall amplitude, which in the pivot frequency is given by |$\alpha _{\rm Tot} = \log _{10} (10^{\alpha _{\rm Fid}} + 10^{\alpha _{\rm PBH}})$|⁠. In particular, to assess the significance of the PBHB contribution, we check whether, for each PBHB population, the value of α_Fid, with its error band, is not compatible at some σ-level (from 1 to 3) with α_Tot, i.e.

where |$\sigma _{\alpha , \rm Tot}$| and |$\sigma _{\alpha , \rm Fid}$| are the FIM errors on α_Tot and α_Fid, respectively, and n ∈ {1, 2, 3}. For reference, we report that for α_Fid = α_SOBHB ≃ −12.02 at f_* = 0.01Hz, we have |$\sigma _{\rm \alpha , Fid}^{4 \rm yr} \simeq 8.44 \times 10^{-3}$| and |$\sigma _{\alpha , \rm Fid}^{10 \rm yr} \simeq 5.34 \times 10^{-3}$| at 68 per cent confidence level after marginalizing over the error on β. Such findings are compatible with those of ref. Babak et al. (2023).

As in the previous section, we conclude by discussing the limitations of our analysis. Analogously to Earth-based detectors and following similar lines of reasoning, we do not include the uncertainties in the population parameters in our analyses.¹⁰ A further crucial approximation is that the FIM formalism gives accurate estimates for the uncertainties in determining the model parameters. This approximation holds in the limit where the log-likelihood for the model parameters is sufficiently Gaussian around |$\vec{\theta }_0$|⁠, which should be quite accurate for the specific injections considered in this work (Babak et al. 2023).

4 RESULTS AND DISCUSSION

We now present the results obtained by applying the proposed methodology to the PBHB population scenario introduced in Section 2. We describe the semi-analytic results achieved by considering the LISA and LIGO |$\rm A^+$| detectors in Section 4.1, while those obtained with LISA and ET in Section 4.2. Given these results, we select a set of benchmark points in the PBHB population parameter space. On these benchmarks, we assess the robustness of the semi-analytic analysis by incorporating spin and sky location parameters, as well as considering realization dependence effects. The results of this assessment are presented in Section 4.3. The overall analysis is conducted on the PBHB parameter space summarized in Table 1, assuming 1 and 10 yr of continuous observations for the Earth-based detectors, and 4 and 10 yr of continuous measurements for LISA.

4.1 Detectability of PBHB populations using LISA and LIGO A⁺

Fig. 2 summarizes the outcome of the semi-analytic analysis showcasing the synergy between LIGO A⁺ and LISA in the slice {μ_LN, ε} of the PBHB parameter space. The third, free parameter of the PBHB model, σ_LN, is set as specified in the title of each panel, which also includes the value of |$\rm T^{\rm LIGO \, A^+}_{\rm Obs}$| adopted in the LIGO A⁺ resolvable-event measurement. For reference, each panel includes the curves corresponding to |$f_{\rm PBH} = 10^{-4}, \, 5 \times 10^{-4}, \, 10^{-3}$| (see Appendix  B) marked as indicated in the legend.¹¹ For each parameter point, we compute the integral in equation (11) and look for the smallest value of z, say |$\bar{z}$|⁠, such that the condition in equation (12) is satisfied. The value of |$\bar{z}$| sets the colour in all these plots. White areas correspond to the parameter regions where the condition in equation (12) is not satisfied at any z. In the colourful region, therefore, LIGO A⁺ would recognize the presence of the PBHB component in the SMBHB population (within the assumptions of the semi-analytic analysis). Moreover, LISA would reach the same conclusion in the PBH parameter regions above the magenta and light brown lines, which are computed by evaluating the condition in equation (23) (see the legend for the LISA observation time, |$\rm T^{\rm LISA}_{\rm Obs}$|⁠, and the sigma level, n, corresponding to each line). In the tiny grey areas appearing in the top right corner of some panels, the overall SGWB signal violates the current LVK upper bound on the SGWB amplitude (Abbott et al. 2021a). Finally, the crosses are the benchmark points that we select for the tests in Section 4.3.

First of all, by scrutinizing Fig. 2, we understand that LIGO |$\rm A^+$| will observe events from the LN PBHB populations only up z ≲ 2. This is both a consequence of the detector’s sensitivity, which only allows for detecting events up to z ≲ 3 (see Section  D), and of the averaging over the spin and angular variables, which slightly lowers the maximal redshift reach compared to, e.g. an optimally located source. Thus, LIGO |$\rm A^+$| will not be able to resolve events in the range where the PBHB population naturally dominates over the SOBHB population, i.e. at redshift higher than the SFR peak. As a consequence, either the fraction of the PBHB population is relatively high at low redshift, or LIGO |$\rm A^+$| will not be able to detect its presence. The additional information provided by the SGWB amplitude measured by LISA might provide an invaluable tool to break the degeneracy among different population models. In particular, the SGWB measurement proves to be quite effective for probing models predicting very narrow peaks at low masses or very broad peaks with small values of ε. On the other hand, for sufficiently large values of μ_LN, populations with narrow mass distributions are more easily detectable with Earth-based detectors. The motivation is that the SNR decreases for increasing mass ratio (q = m₁/m₂). For narrow mass distributions, the two PBHs of each binary are more likely to have similar masses, which, on average, increases the typical event SNR. Moreover, PBHB populations with extremely narrow mass distributions, located at either too small or too large masses, will not be detectable since they generate signals that are either too feeble (the GW amplitude grows with the mass of the binary) or outside the detector’s frequency band (higher masses generally coalesce at lower frequencies). Examples of these effects are visible in, e.g. the left and right top panels of Fig. 2. In particular, in the top left panel, the minimum in the border of the colourful area at |$\mu _{\rm LN} \simeq 70 \, \mathrm{M}_{\odot }$| originates from the interplay between these two effects. Notice also that the shape of the border of the colourful region before such minimum tracks the behaviour of the fiducial mass function, which, as discussed in Appendix  A, is assumed to be a power-law + peak model (Abbott et al. 2023b). Such a shape is expected in the limit σ_LN → 0, and gets smoothed out by increasing σ_LN, as Fig. 2 shows.

Comparing the panels with |$\rm T^{\rm LIGO \, A^+}_{\rm Obs} = 1$| yr with those with |$\rm T^{\rm LIGO \, A^+}_{\rm Obs} = 10$| yr is also interesting. First, we see that, for Earth-based detectors and LISA, increasing the effective observation time generally enhances the detection prospects but not significantly affect the qualitative behaviour of the results. Indeed, for Earth-based detectors, increasing the observation time does not affect the single events detection, but rather, it only increases their number, improving the overall statistics. Similarly, for LISA, increasing the observation time does not impact the overall SGWB amplitude,¹² but rather, the accuracy of its measurement.

Overall, Fig. 2 leads to conclude that, for PBHB mass functions peaked at masses maximizing the SNR, |$\mu _{\rm LN} \simeq 70 \, \mathrm{M}_{\odot }$|⁠, PBH populations with even tiny ε (ε ∼ 10⁻³) are within the reach of LIGO |$\rm A^+$|⁠. While this quickly degrades as PBHB mass functions broaden, this behaviour is not as marked for the prospect of SGWB detection with LISA. As a consequence, parts of the parameter space that can be hardly probed with Earth-based detectors can still be accessed with LISA.

4.2 Detectability of PBHB populations using LISA and ET

We proceed by discussing how the PBHB detection prospects improve when ET replaces LIGO |$\rm A^+$|⁠. Fig. 3 shows the results corresponding to this scenario. It also displays the benchmark points identified in the previous section to highlight their detectability with ET.

As discussed in Section  D, ET can resolve events at z ≳ 10, i.e. well beyond the SFR peak where the SOBHB population quickly drops. As a consequence, ET has way better prospects of identifying PBHB population departures from the SOBHB fiducial model. This fact is manifest in Fig. 3. Moreover, by comparing Fig. 2 and Fig. 3, we see that the detection prospect of ET has less dependence on μ_LN than the one of LIGO |$\rm A^+$|⁠. This effect originates from the improved sensitivity, which, for the mass and redshift ranges considered in this work, leads to less pronounced selection effects in ET, compared to LIGO |$\rm A^+$|⁠. Indeed, all panels of Fig. 3 indicate that selection effects due to the binary mass only affect the very low end of the mass range. We can thus conclude that, as long as the PBHB population will produce a sufficiently large number of events (i.e. larger than the Poissonian 3σ expected for the fiducial population) at high redshift, ET will measure a significant excess in the number of events.

Despite the great increase in the detection ability of ET compared to LIGO A⁺, the SGWB measured by LISA still brings additional information. The main motivation for this claim is that events with very large masses (outside the range of our plots) would merge at too low frequencies to be detected with ET, but would still contribute to the SGWB amplitude. However, a similar argument could also hold for different redshift distributions predicting a few events at low redshift and many more events at very high redshift. Hence, we stress, once again, the synergy between individual events and SGWB for constraining population models.

While, beyond the improvements just underlined, most of the comments explained in Section 4.1 for LIGO |$\rm A^+$| remain valid for ET, we remark that following the methodology introduced in Section 3.1, there are some regions where the PBHB population turns out to be detectable with smaller values of ε in LIGO |$\rm A^+$| than in ET (cf. Fig. 2 and Fig. 3). This might seem counterintuitive, given that ET has better sensitivity. Indeed, this is an artefact of our choice for the detectability criterion, defined in equation (12), and of the selection function for the different GW detectors. With our approach, if the fiducial population produces fewer resolvable events at a given redshift, fewer events are required from the PBHB population to satisfy equation (12). In particular, since LIGO |$\rm A^+$| selects very few events from the fiducial population, it might be easier for a PBHB population with suitable properties (i.e. with a narrow mass function centred at the right value to optimize the SNR at LIGO |$\rm A^+$|⁠) to satisfy equation (12). However, a proper population analysis keeping track of both the redshift and mass distribution (see e.g. ref. Franciolini et al. 2022b) would reveal this feature.

4.3 Analysis of the PBHB population benchmark points

In this section, we perform a more accurate analysis of the resolvable sources with LIGO |$\rm A^+$| and ET focusing on the benchmark points shown in Fig. 2 and Fig. 3. The goal is to quantify the inaccuracies caused by the SNR_avg approximation and the omission of the realization dependencies. The rationale is detailed in Section 3.1.

Let us start by commenting on the choice of our benchmark points. While details are summarized in Table 2, qualitatively these points exhibit the following features:

• Point 1 is chosen to be detectable with LIGO |$\rm A^+$|⁠, but not detectable with LISA/ET, according to the methodology assumed in this work.

• Point 3 leads to signatures in LISA, LIGO |$\rm A^+$|⁠, and ET. Points 2 has the same μ_LN of Point 3 but smaller ε, while Point 4 has the same ε but smaller μ_LN. Both Point 2 and Point 4 are marginally detectable with LISA and LIGO |$\rm A^+$|⁠, but well testable with ET.

The results obtained on these benchmark points using LIGO |$\rm A^+$| and ET are shown in Fig. 4. The three left panels of the figure show the results for LIGO |$\rm A^+$|⁠, while the right panels contain the results for ET. The three rows correspond to three different techniques (with increasing levels of accuracy) for the evaluation of the (expected) number of resolvable sources for a given population. The lines in the top row are computed using the semi-analytic method described in Section 3, and, in particular, by evaluating equation (6). The error bands are estimated from these numbers by assuming a Poisson distribution. Since this procedure is the one used to generate the coloured regions in Fig. 2 and Fig. 3, the population event profiles appearing in these top panels are precisely those establishing the ‘minimum z’ of the benchmarks in those figures (Table 2 quotes such redshifts).

The central panels show the resolvable event distributions determined from the generated catalogues, but still using the SNR_avg approximation. This approximation is the same that we use to get the top-panel results. By comparing the top and central panels, we can assess the impact of realization dependence on our results. As expected, we find consistency in regions with many resolvable sources (small z) and deviations in the low-statistics regime (large z).

Finally, the bottom panels show the results obtained using the catalogues and the complete expression for the SNR, consistently including all waveform variables. Deviations between top (or central) and bottom-panel results manifest at higher redshift. This behaviour is expected since small differences in the source parameters can move borderline sources inside or outside the detection threshold. Thus, including all parameters in the SNR evaluation makes results more dependent on the realization effects in the low-statistics regime. The figure shows that this effect is more pronounced for LIGO |$\rm A^+$|⁠, which has less resolvable sources at high redshift, than for ET. Despite this effect, we observe general good agreement between numerical and semi-analytical results.¹³

The actual number of PBH resolvable events per se is a meaningful quantity. Fig. 5 shows them. Concretely, it displays their redshift evolutions and their Poissonian uncertainties when these quantities are computed through the semi-analytic method. It also includes the uncertainty |$\Delta ^{\rm Res, Fid}_z$| (at 1 and 3σ-level) for the fiducial population. The figure gives an insight into how the realization dependencies on the PBH populations would influence the conclusions reached via our semi-analytic analysis. We can see that the impact is non-negligible but does not change our main conclusions.

Before concluding, we further scrutinize the relative differences due to the SNR approximation and the realization dependencies. For this purpose, using the same notation introduced in Section 3.1, we define the following quantities:

where ‘hat’ quantities are computed using sky and spin averages, and ‘bar’ quantities use the sky and spin information. Fig. 6 shows their values as a function of redshift for the fiducial population in blue and red, respectively. We see that before reaching the low-statistics regime, the blue curves are always smaller than |$10{{\ \rm per\ cent}}$| for both LIGO A⁺ and ET. This proves a quite good agreement between the semi-analytical approach and the generated catalogues. On the other hand, the red curves reach higher values (⁠|$\simeq 30 / 40{{\ \rm per\ cent}}$|⁠), so the SNR approximation can induce larger errors than the realization dependence. Nevertheless, both the SNR approximation and the realization dependence are sufficiently small not to jeopardize the main conclusions that one would reach via the semi-analytic method.

5 CONCLUSIONS

In this paper, we have discussed the prospects of detecting potential PBHB populations with future GW detectors. For this purpose, we have assumed a fiducial population in agreement with the GWTC-3 results and added PBHB populations with different merger rates and mass distributions to test whether these would lead to observable signatures in LIGO |$\rm A^+$|⁠, ET, or LISA. Using Earth-based detectors, we have checked how the number of resolvable events changes in the presence of PBHB populations. In particular, we have evaluated this analytically and tested our results by simulating event catalogues to assess the impact of low statistics on the analytic results. We have generally found good agreement between our semi-analytical and numerical results. Beyond that, we have evaluated the increase in the amplitude of the SGWB arising from the PBH contribution and tested its detectability with LISA using a FIM approach. We found that the information LISA will bring might be significant to test whether the SGWB is due to SOBHBs only.

For all models considered in this work, we have found sizable regions of the parameter spaces where the PBHB populations will lead to significant variations in the number of detectable events in LIGO |$\rm A^+$| and ET (with ET performing better in most cases) with respect to SOBHB expectations. However, it is worth stressing that detecting events at high redshift does not imply that it will always be possible to infer their redshift accurately (Ng et al. 2022a; Mancarella et al. 2023). Moreover, we have found that for all models considered in this work, there are sizable parts of the parameter spaces leading to an increase in the SGWB amplitude that would be detectable with LISA. Interestingly, since the SGWB integrates over all masses, the SGWB measurement can also test populations with very low and very large masses, generating signals beyond the reach of future Earth-based detectors. Indeed, different GW detectors probe complementarily distinct parts of the parameter space. In particular, our results highlight three different regimes:

• Signatures in Earth-based detectors and LISA: This can happen if the PBH population becomes abundant (but still statistically marginal in present LVK observations) around (or after) the SFR peak so that the number of individual events does not decline at z ≃ 2. Simultaneously, the SGWB at LISA exceeds the SOBHB prediction.

• Signatures in Earth-based detectors only: The PBH population grows very slowly with z and becomes sizable only at high redshifts. In this case, the signals from unresolved sources are faint, and their contribution to the SGWB at LISA is not sufficiently strong to modify the SOBHB prediction significantly. While deviations in the merger rate might be appreciated, their statistical significance might not be sufficient to pin down the presence of a secondary population unless some additional features are found in, e.g. the mass distribution of the observed population.

• Signatures in LISA only: The PBH population has very small or very large masses, and the signals are not detectable with Earth-based detectors. LISA observes an SGWB amplitude incompatible with the value predicted using the SOBHB population measured by Earth-based detectors.

Overall, we conclude that the considered measurements from Earth-based detectors and LISA will generally be complementary, and our understanding of the BH population that we observe in our Universe will improve if we use their synergy. This conclusion comes with no surprise, knowing the underlying differences between the properties of the signals these detectors will probe. Furthermore, we generally observe that the dependence of the LISA SGWB on the population parameters scales differently than the distribution of resolvable sources that Earth-based interferometers will detect. This fact implies that, in general, SGWB measurements will help Earth-based detectors improve the constraints on the BH population parameters that we observe in our Universe.

As discussed in Section 3, our analysis does not include errors on the fiducial population parameter, which are currently quite broad on some of the most influencing parameters and would impact our results significantly if extrapolated to the volume that LIGO |$\rm A^+$| and ET detector will probe. However, the open codes presented in our GitHub repository (Marcoccia 2023a) can be readily updated when the new results of future inference papers, e.g. by the LVK collaboration, come out. With improvements in the LVK network, we expect more (and more accurate) detections, which will reduce the uncertainties on the population parameters, making the analysis much more reliable.

This study assumes that the two populations do not interact with one another, i.e. mergers only involve BHs drawn from the same population. This assumption impacts both the number and the properties of the mergers. By dropping this assumption and assuming the two populations to have sufficiently different mass ranges, it would be possible to enhance the number of Extreme Mass Ratio Inspirals in the LISA band significantly (Babak et al. 2017; Gair et al. 2017; Guo & Miller 2022; Mazzolari et al. 2022).¹⁴ Thus, determining the abundance of these objects can also be used to further constrain the eventual presence of PBHB populations on GW detectors with LISA-like frequency range. Moreover, keeping track of the number of resolvable sources in different frequency bands (e.g. BHs with masses higher than the ones considered in this work would merge in the LISA frequency band) could also provide a possible tracer for the presence of PBHB population. Finally, including the SGWB measurement at the Earth-based detectors could improve the results presented in this study for two main reasons. First, detecting the SGWB at different frequency bands will decrease the uncertainties on its amplitude, hence improving the chances of detecting deviations from the expected value. Secondly, if the PBH mass function is narrow and peaked in the stellar-mass range, the PBH could modify the SGWB shape introducing deviations from the standard power-law behaviour in the LVK/ET frequency range, which could be used as a further constraint on the PBH sub-population properties (Chen et al. 2019; Kapadia et al. 2020; Bavera et al. 2022; Franciolini & Pani 2023).

We conclude by commenting on alternative methods to assess the detectability of PBHB populations beyond the ones considered in this work for both Earth-based and Space-based detectors. It would be possible, in principle, to adopt other criteria similar to the one we have introduced. For example, since we expect PBHBs to become relevant at high redshift, variations on the cumulative number of resolvable sources predicted after a given redshift could provide a viable alternative. While maintaining less information on the source distribution, this approach would be less sensitive to the error in the inference of the source distance (Ng et al. 2022a; Mancarella et al. 2023). Finally, a fully consistent population analysis based on hierarchical Bayesian modelling (see e.g. Wong et al. 2021; Chen, Yuan & Huang 2022a; Franciolini et al. 2022b) would be computationally more expensive but, keeping its subtleties under control, more accurate than the criteria discussed here. Thus, we deem it worth exploring these (and possibly other) methodologies in future works.

ACKNOWLEDGEMENTS

We thank Gabriele Franciolini for several useful comments on a draft of this work. MP acknowledges the hospitality of Imperial College London, which provided office space during some parts of this project.

Funding. PM acknowledges the Cost action CA16104 for financing the STSM to London that led to the start of this work. GN is partly supported by the grant Project. No. 302640 funded by the Research Council of Norway.

DATA AVAILABILITY

The authors confirm that the data supporting the findings of this study may be generated using the open codes presented in our GitHub repository (Marcoccia 2023a). The catalogues for both the fiducial and PBH populations may be generated using the open code presented in our Github repository (Marcoccia 2023b), which follows the procedure described in ref. Babak et al. (2023). The raw catalogue realizations used to produce the results of this paper are available from the corresponding author, upon reasonable request.

Footnotes

REFERENCES

APPENDIX A: THE SOBHB FIDUCIAL POPULATION

In this appendix, we detail the fiducial SOBHB population model we adopt throughout the analysis. In most aspects, we follow the approach of ref. (Babak et al. 2023) and rely on the master equation in equation (1). We proceed by clarifying the functional forms of the quantities appearing in such an equation.

As already discussed in the main text, the current LVK data put tight bounds on the SMBHB population properties up to z ∼ 0.5 (Abbott et al. 2023a). Few events have been detected at redshift 0.5 ≲ z ≲ 1, but they are too rare and/or poorly reconstructed to impose strong constraints (Abbott et al. 2020; O’Brien et al. 2021). Despite these caveats, current data are compatible with a population of SOBHBs with a merger rate behaving as

for z ≲ 0.5, with |$R(z = 0.2) = 28.3^{+13.9}_{- 9.1} {\rm Gpc}^{-3} {\rm yr}^{-1}$| and |$\kappa = 2.9^{+ 1.7}_{- 1.8}$| (Abbott et al. 2021b, 2023a). At higher redshift, R(z) has to keep track of the stellar-formation origin of the binaries and, to some degree, resemble the SFR. Thus, consistently with Babak et al. (2023), we choose the Madau-Dickinson phenomenological profile (Madau & Dickinson 2014; Mangiagli et al. 2019) with a negligible time delay.¹⁵ Such a choice leads to

where R₀ is set so that R(z = 0.2) matches the measured value.

For what concerns the mass distribution, we adopt the power law + peak scenario (Abbott et al. 2021b, 2023a)

with C_mass being a normalization constant and q being the mass ratio q = m₂/m₁. The functions π₁ and π₂ read as

and

where C_q(m₁) is a normalization function, |$\mathfrak {S}$| is a smoothing function for the low mass cutoff, and |$\mathfrak {P}$| and |$\mathfrak {G}$| are respectively a normalized power law and a normalized Gaussian distribution

with α being the spectral index of the power law, μ_m and σ_m being the mean and width of the Gaussian, and C_PL and C_m being normalization. The smoothing function |$\mathfrak {S}$| imposes a smooth cutoff for low masses, rising from 0 to 1 in the interval [m_m, m_m + δ_m]

with

so that, by construction, we have m ≥ m_m. The high end of the mass range does not have an explicit cutoff, but large masses are statistically suppressed. In practice, we set m_M = 100M_⊙, which is slightly higher than the values used in the LVK analysis (Abbott et al. 2023b), to take into account possible higher mass events of astrophysical origin.¹⁶All the hyperparameters entering the mass distribution equation (A3) are fixed at the central values of the LVK analysis outcome reported in Table A1.

The spin distribution is a product of two different PDFs, one for the spin amplitudes and one for the spin tilts. The former reads as (Abbott et al. 2021b, 2023a):

where B(α_a, β_a) is a Beta function that guarantees the appropriate normalization of the PDF. The α_a and β_a are positive constants defined through

where E[a] and Var[a] are set in Table A2. We stress that the spin amplitudes of the two black holes are independent of one another. On the other hand, the PDF spin tilt distribution reads

which is a mixture between an isotropic and a truncated Gaussian distribution centred in cos (t_i) ≈ 1. The σ_i, and ζ parameters are specified in Table A2.

Finally, for cosmology we assume the ΛCDM model where the Hubble parameter is , with h = 0.678 being its dimension-less|$H_0 = h \times 100 \, {\rm km\,(s \, Mpc)^{-1} }$| value, and Ω_m = 0.3 and |$\Omega _\Lambda = 0.7$|⁠.

APPENDIX B: PBH CONTRIBUTION TO THE DARK MATTER RELIC ABUNDANCE

While PBHs behave as cold DM and could, at least in principle constitute a sizable amount of the presently observed DM, their abundance in the stellar mass range is tightly constrained (Nakamura et al. 1997; Ioka et al. 1998; Tisserand et al. 2007; Ricotti et al. 2008; Wyrzykowski et al. 2011; Ali-Haïmoud et al. 2017; García-Bellido 2017; Villanueva-Domingo et al. 2021). We define f_PBH ≡ Ω_PBH/Ω_DM, the ratio between today’s PBH and DM energy densities in the Universe. For a given PBH population model, the parameters ε and f_PBH can be explicitly related. We obtain their relationship in the case that our PBH population, which we derive from phenomenological models, is dominated by binaries that gravitationally decoupled from the Hubble flow before the matter-radiation equality (Hütsi et al. 2021; Bagui et al. 2023; Franciolini et al. 2023).

In equation (5), we modulate the PBHB merger rate R_PBHB(z) through the parameter ε. An alternative way to write R_PBHB(z) is (Hütsi et al. 2021; Franciolini et al. 2023)

where η = m₁m₂/(m₁ + m₂)² and |$\mathcal {S} = \Phi _{\rm LN}(m_1) \, \Phi _{\rm LN}(m_2) S$|⁠. The function Φ_LN is given in equation (7). The function S is a suppression factor accounting for environmental effects that slow down the binary formation or favour their disruption. It can be approximated as¹⁷

with

Here σ_M ≈ 0.004 represents the rescaled variance of matter density perturbations at the time the binaries form, Γ denotes the Euler gamma function, U(a, b, z) is the confluent hypergeometric function, while <m > and <m² > are the first and second momenta of the PBH mass PDF.

The comparison of equation (5) to equation (B1) yields

while Fig. B1 shows the values of f_PBH that arise in the parameter regions of the PBH models considered in our analysis.

APPENDIX C: IMPACT OF THE FACTORIZATION IN THE MASS FUNCTION DISTRIBUTION

In the PBH literature, the PBHB mass PDF is typically expressed as

while in our analysis, we have approximated it as

The motivations for this approximation are both to improve the computational time of the analysis, as well as to give a clearer qualitative description of the properties of the PBH populations that we are considering. The main effect of the terms we neglected is to drag the PBHB total mass toward lower values compared to the cases considered in this work. However, we expect these corrections to be relevant only for populations with high σ_LN, while for populations with small σ_LN the exponential decay of the log-normal PDF dominates the generated masses distribution. In this appendix, we quantify the impact of neglecting such mixing terms in our analysis, and in particular, for the results in Fig. 2 and Fig. 3.

Concerning the SGWB analysis, we have computed numerically the PBHB SGWB amplitude with the exact and approximated mass PDF expressions in the scenario with σ_LN = 1. We find that the mixing term lowers the SGWB amplitude by less than |$50{{\ \rm per\ cent}}$| even for high values of μ_LN. This implies that the SGWB amplitude remains within the same order of magnitude, with a minimal variation of the SGWB lines presented in Figs 2 and 3.

To assess the impact of our approximation on the results for the resolvable events, we show, in Fig. C1, the results on the benchmarks Point 1 and Point 3 (the former exhibits a small value of σ_LN, the latter exhibits a maximal value of σ_LN and high number of events) without neglecting the terms in equation (C1). In the panels, we label as |$\tilde{1}$| and |$\tilde{3}$| the results including the mixing terms. Fig. C1 shows that the correction, which tends to decrease the SNR of the events, is not noticeable for the ET detector. This comes with no surprise, as in ET the PBHB populations tend to have an SNR much higher than the assumed detection threshold in the considered redshift interval. Even for LIGO A⁺, the mixing terms do not affect significantly the results for benchmark Point 1, which has a narrow mass distribution. On the other hand, for Point 3, which corresponds to a broad mass function, the dampening in the expected number of resolvable sources is appreciable but only in the high-redshift region of the plot, in particular much above the redshift at which the PBHB population emerges from the SOBHB fiducial model (cf. Fig. 2 and Fig. 3).

Finally, in Fig. C2 we compare the impact of this approximation, for the case of Point 3, with the other approximations that we considered in Fig. 4. It results that, even for σ_LN = 1, this effect is subdominant when compared to the impact of realization and SNR approximation errors.

APPENDIX D: DETECTOR CHARACTERISTICS

In this study, we consider LIGO A⁺ and ET as representatives of upcoming and future Earth-based interferometers and LISA as a reference for the first generation of space-based GW detectors. The precise timeline and operational durations of these instruments are uncertain. Nevertheless, it is reasonable to anticipate that LIGO A⁺ will operate for several years before the early/mid-2030s when ET and LISA are expected to commence to acquire data. In the lack of a well-defined progress plan, we consider a couple of somewhat extreme timeline scenarios, believing that the actual future will likely fall somewhere in between. Concretely, we analyse 1 and 10 yr of data for LIGO A⁺ and ET, and 4 and 10 for LISA. We leave it to the knowledge of the future reader to estimate which scenario the future will tend to and in which order each detector and its measurements will arrive.

D1 Earth-based interferometers

The location of the two LIGO A⁺ detectors are set to be in the Livingston (N 30°330′, W 90°460′) and Hanford (N 46°270′, W 119°240′) sites. Regarding the sensitivity, we use the curve described in ref. (LIGO 2018), with frequency range [5, 5000] Hz. For ET, we assume the location proposed in the Sos Enattos mine in the Lula area (N 40°260′, E 9°260′) with the ET-D-sum sensitivity in the frequency range [0.1, 10⁴] Hz (Hild et al. 2011; ET 2018). However, our resolvable event analysis is nearly independent of the precise detector sites. The expected horizon distance for these detector configurations, w.r.t the BHB populations considered in this paper, is presented in Fig. D1

D2 LISA

LISA will be the first interferometer in space. The detector will consist of three satellites orbiting around the Lagrange point L5. For our analysis, we assume mission adoption sensitivity in the frequency range [3 × 10⁻⁵, 0.5] Hz (LISA 2018). In the following, we describe the 2-parameters instrument noise model (Amaro-Seoane et al. 2017; Babak, Petiteau & Hewitson 2021) based on the results of the LISA Pathfinder mission, as well as the latest laboratory test. In particular, we report the LISA sensitivity in the Time Delay Interferometry (TDI) channels A and T.¹⁸ For more details on the noise model and the TDI construction, see e.g. ref. (Flauger et al. 2021; Hartwig et al. 2023).

The noise in LISA is a combination of two main components: Test Mass (TM) acceleration noise and Optical Metrology System (OMS) noise. The power spectra P_TM, P_OMS, for these two components are

where c is the light speed and the two noise parameters A and P control the amplitudes of the TM and OMS components, respectively. The total noise spectral densities in the TDI A and T channels read

with |$L = 2.5 \, \times 10^{9}$|m is the LISA armlength.

Given the noise power spectra, the strain sensitivity (for a generical channel ij) is defined as

where R_ij(f) is the (quadratic) response function, mapping incoming GW signals onto the TDI data stream. The response can be further expanded as a purely geometrical factor, |$\tilde{R}_{ij}(f)$|⁠, times TDI-dependent terms. While |$\tilde{R}_{ij}(f)$| should be evaluated, approximate expressions for the A and T channels read

Since |$\tilde{R}_{\rm TT}$| is strongly suppressed at low frequencies with respect to |$\tilde{R}_{\rm AA}$|⁠, the T channel is typically assumed to be signal insensitive. It is customary to express the noise in Ω units using

In the analyses presented in this work, we assume the fiducial values for the noise parameters to be A = 3, P = 15 with 20 |${{\ \rm per\ cent}}$| Gaussian priors.

APPENDIX E: ANALYTICAL DERIVATION OF THE SGWB FROM A POPULATION OF MERGING OBJECTS

In this appendix, we summarize the analytical derivation of the SGWB sourced by the PBHB and SOBHB populations. Further details on the derivation can be found in refs. (Phinney 2001; Sesana, Vecchio & Colacino 2008) for details. We focus on the case of the LISA detector, for which we know that the number of resolvable sources is small enough (∼10 sources for |$\rm SNR_{\rm tresh}= 8$|⁠) to prevent dampenings of the SGWB signal in the high-frequency region of the LISA sensitivity band (Babak et al. 2023).

Following the approach of refs. (Hawking & Israel 1989; Phinney 2001), the characteristic spectral strain of the SGWB can be defined as

Here dn/dz and dE_GW/df_r are, respectively, the comoving number density of the sources, and the (redshifted) energy per log frequency interval produced by each source. The latter, in the circular-orbit approximation, reads as

with |$\mathcal {M}$| being the chirp mass of the source. The former is a function of the statistical properties of the source population. It can be written as

where R(z) is the merger rate of the population, |$P[\mathcal {M}(m_1, m_2)]$| is the probability of having a chirp mass |$\mathcal {M}$| from the m₁ and m₂ mass PDFs, and dt_r/dz is a drift term defined as

By putting all these elements together, one finds

This can be rephrased in Ω units as