© The Authors 2024

1. Introduction

The Laser Interferometer Gravitational-wave Observatory (LIGO) (Aasi et al. 2015), Virgo (Acernese et al. 2015), and KAGRA (Akutsu et al. 2019) detectors are progressively uncovering the features of the population of merging stellar-mass binary black holes in our Universe (Abbott et al. 2021a, 2023a). As observing runs become more and more sensitive, the detection horizon increases and a higher number of gravitational-wave (GW) events are positively identified as binary mergers. The events observed so far lie at distances z ≲ 1 (Abbott et al. 2023b), while the vast majority of binary mergers are expected to lie well beyond this horizon, as suggested by both theoretical expectations for the merger rate redshift evolution (Dominik et al. 2015; Belczynski et al. 2016; Rodriguez et al. 2018) and its inferred trend from GW data through the Third Gravitational-Wave Transient Catalogue (GWTC-3) Abbott et al. (2023a). This knowledge has motivated studies and the inference of the sub-threshold unresolved collection of binaries, treated as an overall gravitational-wave background (GWB) signal.

Works forecasting the GWB from compact binaries populate the literature. Before the first GW detections, these were based on theoretical models of the binary population (Phinney 2001; Regimbau 2007, 2011), while more recently GW data-informed projections (Abbott et al. 2021c) have become a benchmark for the GW community. An important distinction to make is between estimates of a specific realisation of the background, for example used in mock-data challenges (Meacher et al. 2014), and estimates of the ensemble average of the background, which corresponds to the expectation value of the background amplitude targeted by GWB searches (Zhu et al. 2011; Wu et al. 2012). In both cases, the calculation involves the (expected) distributions for the individual binary parameters, such as mass, distance, and event rate distributions. These population models are often taken to be simple parametric functions with a well-defined set of hyperparameters, which are fixed to fiducial (or assumed) values for the background calculation. Recalculating the background signal for varying population hyperparameters can become computationally intensive when employing large sample sets.

Several applications in the literature require marginalising the background signal over possible population configurations, including forecasting studies, such as those presented by the LIGO-Virgo-KAGRA (LVK) collaboration in Abbott et al. (2021c) and Abbott et al. (2023a), and inference analyses in Callister et al. (2020), Abbott et al. (2021c), and Turbang et al. (2024), among others. As the interest in this type of work grows, there is a need for efficient and flexible background estimation procedures. In this paper we present a method to efficiently carry out these calculations, called popstock, that we have made available to the community as an open-access code base. The novelty of our approach is in the design of a reweighting technique that allows us to sample the binary parameter probability distributions and evaluate a corresponding set of waveform approximants only once, enabling an extremely efficient re-estimation of the GWB when varying population hyperparameters. An analogous reweighting approach was previously implemented in Turbang et al. (2024). Furthermore, we implement the use of waveform templates imported from the LIGO Scientific Collaboration Algorithm Library (LAL) LIGO Scientific Collaboration (2020). To improve efficiency, previous codes have employed analytic waveform approximations directly embedded in the codebase; with our CPU-optimised spectral calculations and our GPU-optimised reweighting technique, we are able to support a much broader range of waveforms through the commonly used Python library for GW analysis, Bilby (Ashton et al. 2019).

We used popstock to investigate the impact of the uncertainty on the binary population on the detection prospects of the GWB with ground-based interferometers. In particular, we made a comparison with and extended the work done in Abbott et al. (2023a) to include more uncertainties on the redshift distribution of sources. We also employed the package to probe the effect of waveform choice on the estimation, and whether this is entirely degenerate with the population uncertainty. We find that the expected background amplitude can be significantly boosted when admitting higher mass mergers and higher rates of mergers at low redshift, within current population uncertainties. We further find that the choice of waveform can have noteworthy effects on background signal estimates; however, these do not dominate the current population uncertainties.

This paper is organised as follows. In Sect. 2 we introduce the theoretical aspects behind compact binary GWB calculations. In Sect. 3 we define the analytic models used here to describe the distributions of black hole masses and distances. In Sect. 4 we introduce popstock, our new package for GWB calculations, outlining its key functionalities. In Sect. 5 we present our background calculations and investigations, and finally we summarise our conclusions in Sect. 6.

2. The compact binary GW background

The amplitude of the GWB signal is parametrised by the fractional energy density spectrum emitted by GWs throughout the Universe, Ω_GW(f) (Phinney 2001),

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{GW}}(f)=\frac{1}{{\rho }_{\mathrm{c}}}\frac{\mathrm{d}{\rho }_{g}w}{\mathrm{d}lnf},\end{array}$$(1)

which is normalised by the critical energy density of the Universe ρ_c = 3c²H₀²/8πG. This is in general the total energy density contributed by GWs throughout the Universe, and is not restricted to sub-threshold signals. Calculating a residual background requires the definition of a detection cutoff, which is detector-dependent¹, while here and elsewhere Ω_GW is considered to be an astrophysical property of the Universe.

While the GWB spectrum is detector-independent, it is useful to employ the detector frame to perform calculations. This allows us to immediately relate the intrinsic signal amplitude with the measured signal in GW detectors. As shown for example in Meacher et al. (2015), and Regimbau (2022), the energy density spectrum Ω_GW(f) may be estimated as the average fractional GW energy density present in a detector during an observation time T_obs,

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{GW}}(f)=\frac{f^3}{T_{\mathrm{obs}}}\frac{4{\pi }^2}{3H_0^2}{\displaystyle \sum _i^{N_i}}{\mathcal{P}}_{\mathrm{d}}({\mathrm{\Theta }}_i;f),\end{array}$$(2)

assuming a finite number of GWs received at the detector, N_i. Here 𝒫_d(Θ_i, f) is the Fourier domain unpolarised power in the detector frame (hence the subscript d) associated with a GW with parameters Θ_i, defined as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{d}}({\mathrm{\Theta }}_i;f)={\stackrel{\sim }{h}}_{+}^2({\mathrm{\Theta }}_i;f)+{\stackrel{\sim }{h}}_{\times ,i}^2({\mathrm{\Theta }}_i;f),\end{array}$$(3)

where ${\stackrel{\sim }{h}}_A({\mathrm{\Theta }}_i;f)$ is the Fourier transform of the GW waveform evaluated at Θ_i ². So defined, 𝒫_d has units s⁻².

In the limit of infinite observation time and infinite events, the Ω_GW spectrum approaches its ensemble average, ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$. While the Ω_GW measured by an experiment depends on the specific realisation observed throughout the experiment observation time, ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ only depends on the distributions that describe the GW parameters Θ, which are in turn parametrised via population hyperparameters Λ (Phinney 2001). The ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ spectrum is thus a unique population signature, and is targeted in observations in practice by measuring Ω_GW for large T_obs. The equivalence is easily seen in the limit of large GW numbers as

$$\begin{array}{c}\hfill \frac{1}{T_{\mathrm{obs}}}{\displaystyle \sum _i^{N_i}}{\mathcal{P}}_{\mathrm{d}}({\mathrm{\Theta }}_i;f){\displaystyle \underset{\begin{array}{c}{N}_i\to \infty \\ T_{\mathrm{obs}}\to \infty \end{array}}{\approx }}\frac{\mathrm{d}N}{\mathrm{d}t}{\displaystyle \int \mathrm{d}\mathrm{\Theta }}{p}_{\mathrm{d}}(\mathrm{\Theta }|\mathrm{\Lambda }){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta };f),\end{array}$$(4)

where p_d(Θ|Λ) are the (normalised) detector frame probability distributions for the GW parameters Θ, such that

$$\begin{array}{c}\hfill {\overline{\mathrm{\Omega }}}_{\mathrm{GW}}(\mathrm{\Lambda };f)=f^3\frac{4{\pi }^2}{3H_0^2}R{\displaystyle \int \mathrm{d}\mathrm{\Theta }}{p}_{\mathrm{d}}(\mathrm{\Theta }|\mathrm{\Lambda }){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta };f).\end{array}$$(5)

Here $R\equiv \frac{\mathrm{d}N}{\mathrm{d}t}$ is the total rate of events per unit detector-frame time. It can be convenient to isolate the redshift integral in Eq. (5), assuming redshift is independent from other parameters, and to incorporate the rate in the redshift evolution probability p(z), defining the event rate per unit detector-frame time per redshift shell:

$$\begin{array}{c}\hfill R(z)=Rp(z).\end{array}$$(6)

The ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ spectrum can then be interpreted as an integral over redshift shells of the average GW power present in each shell, analogously to Eqs. (4) and (5) in Callister et al. (2020),

$$\begin{array}{c}\hfill {\overline{\mathrm{\Omega }}}_{\mathrm{GW}}(\mathrm{\Lambda };f)=f^3\frac{4{\pi }^2}{3H_0^2}{\displaystyle \int \mathrm{d}}zR(z|{\mathrm{\Lambda }}_z)⟨{\mathcal{P}}_{\mathrm{d}}(z,\theta ;f)⟩,\end{array}$$(7)

where the ⟨.⟩ brackets imply the GW spectral power samples 𝒫_d are averaged over the ensemble described by the parameter probability distributions, in each redshift shell. To see how the GW spectral power is related to the energy spectrum, for each binary (see Appendix A for a pedagogical derivation).

3. BBH population models

We illustrate the effect of the population model on Ω_GW by adopting two phenomenological mass distribution models used in Abbott et al. (2021a, 2023a). The POWERLAW+PEAK model (PLPP), first introduced in Talbot & Thrane (2018), has been widely adopted in the literature as an astrophysically motivated mass distribution model. The PLPP model consists simply of a truncated power law, motivated by the shape of the stellar initial mass function, and a Gaussian bump (or peak), originally intended to account for a possible overdensity of black holes around a certain mass, as motivated for example by pulsational pair instability effects (Talbot & Thrane 2018; Abbott et al. 2021a, 2023a)³. We list the parameters of the PLPP model in Table 1. In addition to PLPP, we also consider a simpler mass distribution, consisting of a truncated power law with a break at a particular mass. While Abbott et al. (2023a) finds that this broken power-law (BPL) model is disfavoured with respect to the PLPP model, we include it for illustration purposes. The parameters used for the BPL model are described in Table 2.

We also adopt the redshift distribution models commonly used in the literature. For example, Fishbach et al. (2018) introduced a broken power-law model to describe the merger rate as a function of redshift. This model is motivated by the observed star formation rate (SFR) across redshift: a rate rising to and peaking at some redshift (z_peak; Mandel & de Mink 2016; Fishbach et al. 2018) then decaying down to high redshift, where the SFR was much lower. A SFR model commonly adopted in the literature is the broken power-law fit from Madau & Dickinson (2014), which is parametrised in terms of a low-redshift power-law index γ, high-redshift (negative) index κ, and a peak or turnover redshift parameter z_peak:

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{MD}}(z)\propto \frac{{(1+z)}^{\gamma }}{1+{\left(\frac{1+z}{1+z_{\mathrm{peak}}}\right)}^{\kappa }}.\end{array}$$(8)

In the following, as in Abbott et al. (2023a) and Abbott et al. (2021c), we use ℛ_MD to describe the number of mergers per unit comoving volume V_c per unit source-frame time t_s. For brevity, we refer to this model as MD. The parameters employed throughout for the MD model and their default values are shown in Table 3. In reality, the merger rate of BBHs from stellar collapse is a function of the SFR and the delay time between star formation and merger of the remnant Fishbach & Kalogera (2021), Callister et al. (2020), Turbang et al. (2024); this may make the number of mergers per unit comoving volume and source-frame time R(z) deviate from a SFR-like broken power law.

For illustration, we also adopt a less realistic model in which the merger rate is constant across cosmic time, ℛ(z)=R₀. This uniform-in-comoving-volume (UICV) model assumes that the merger rate as measured in the source frame of the emitter, is constant across redshift.

The redshift of individual compact binary merger events in the detector frame is drawn from a probability distribution p(z), which takes into account the comoving volume per unit redshift gradient dV_c/dz, and the redshifting of the rate from source-frame to detector-frame,

$$\begin{array}{c}\hfill p(z)=\frac{1}{1+z}\frac{\mathrm{d}{V}_{\mathrm{c}}}{\mathrm{d}z}\mathcal{R}(z).\end{array}$$(9)

The p(z) functions for the MD and UICM rate evolution are illustrated in the left panel of Fig. 3.

In this paper, unless otherwise specified, we draw uncertainties from the LVK collaboration GWTC-3 population posteriors, published in the data release LIGO Scientific Collaboration, Virgo Collaboration, KAGRA Collaboration (2021) which accompanied the collaboration results Abbott et al. (2023a). The release includes samples from the posterior of population hyperparameters inferred through GWTC-3 (i.e. the population parameters governing the shapes of the mass, spin, and redshift distributions). We use these hyperparameter samples for the corresponding redshift and mass models described above in the analysis that follows. We use the PLPP mass model and MD redshift model as fiducial models for our studies. The PLPP model is considered a good parametric description of the mass spectrum of GWTC-3 Abbott et al. (2023a), also confirmed by non-parametric approaches Callister & Farr (2024), and has been widely used since in the context of GWB estimation (e.g. Zhou et al. 2023; Bellie et al. 2024; Sah & Mukherjee 2023), while the MD model is a well-motivated astrophysical model (Fishbach et al. 2018) and has been already employed in stochastic inference analyses (Abbott et al. 2021c).

4. popstock

We present popstock⁴, a Python-based open-source package for the rapid computation of background spectra such as Ω_GW, for a given realisation of events, and ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$, for a given set of hyperparameters Λ. Other than the standard Python scientific libraries numpy (Harris et al. 2020) and scipy (Virtanen et al. 2020), the main dependencies of the popstock package are astropy (Astropy Collaboration 2022), a core Python library used by astronomers; bilby (Ashton et al. 2019), a very popular Bayesian inference library for GW astronomy; and gwpopulation (Talbot et al. 2019), a hierarchical Bayesian inference package containing a collection of parametric binary black hole mass, redshift, and spin population models.

The popstock package relies on multiprocessing (included in most Python distributions) to parallelise the computation of Ω_GW for large N_i. The GW waveforms required to compute 𝒫_d(Θ_i; f) in Eq. (3) are evaluated at each Θ_i using the Bilby library, which in turn imports LAL LIGO Scientific Collaboration (2020). This allows us to employ a vast array of modern waveforms in our computations.

To compute ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ for a given set of population hyperparameters Λ and a given collection of population models, we directly sample the probability distributions p_d(Θ|Λ) and evaluate Eq. (5) via a Monte Carlo simulation. The accuracy of this evaluation depends on the number of samples employed, as discussed below. This approach is limited by the long sampling and evaluation times of the GW waveforms, and is not an optimal tool for performing in-depth studies of the impact of population uncertainties on the ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ spectrum. Hence, popstock includes a reweighting technique to compute ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ for a new set of Λ parameters without re-evaluating Eq. (5). In the rest of this section, we describe the popstock reweighting technique and probe its efficiency and accuracy.

The popstock repository includes tutorials with usage instructions and simple examples. The package’s performance strongly depends on the computing setup employed (e.g. CPU and GPU availability). We refer to the popstock documentation⁵ for further details.

Fig. 1. Impact of sample variance and reweighting on the ΩGW spectrum. Left: 95% confidence on the spectrum calculated using 105 (106) samples in yellow (green) drawn from a fixed hyperparameter distribution, compared to the 95% confidence on the spectrum including uncertainty on the local merger rate parameter in brown. Unsurprisingly, the sample variance varies significantly with number of samples and when using large sample numbers is subdominant compared to population parameter uncertainty. Right: sample variance from the left panel compared to reweighted estimates of the ΩGW spectrum. The reweighted spectra lie neatly within the sample variance uncertainty bounds, implying that the reweighted spectrum is indistinguishable from a regularly sampled spectrum, with these sample numbers.

4.1. Reweighting methodology

We lay out a simple method to efficiently calculate ${\overline{\mathrm{\Omega }}}_{\mathrm{GW}}$ for different sets of hyperparameters Λ_i describing the (detector frame) target population distributions, p_d(Θ|Λ_i). The integral of (5) above allows the implementation of an importance sampling approach or reweighting,⁶ whereby

$$\begin{array}{cc}\hfill {\displaystyle \int \mathrm{d}\mathrm{\Theta }}{p}_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_1){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta })& ={\displaystyle \int \mathrm{d}\mathrm{\Theta }}\frac{p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_1)}{p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_0)}{p}_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_0){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta })\hfill \\ \hfill & \equiv {\displaystyle \int \mathrm{d}\mathrm{\Theta }}{w}_{0\to 1}(\mathrm{\Theta })p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_0){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta }),\hfill \end{array}$$(10)

where p_d(Θ|Λ₀) is chosen as the fiducial distribution, and w_0 → 1 is the weight required to ‘transform’ between the fiducial distribution and the target distribution, relative to Λ₁:

$$\begin{array}{c}\hfill w_{0\to 1}=\frac{p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_1)}{p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_0)}.\end{array}$$(11)

In practice, this reweighting approach is more efficient than direct Monte Carlo integration when p_d(Θ|Λ₁) is hard to sample from but easy to evaluate. Therefore, we first draw a large set of samples θ from the fiducial population Λ₀ and compute the probability of drawing those samples p_d(Θ|Λ₀). This is stored as the denominator in the weights w. Each time the integral of Eq. (10) is evaluated for some different population Λ₁, it is only necessary to evaluate the probability of those fiducial samples under the target population p_d(Θ|Λ₁), as w is the only term that depends on Λ₁. See also Appendix D of Turbang et al. (2024) for an analogous reweighting approach for the background spectrum calculation. The authors’ approach is technically equivalent; however, they choose to reweight the approximated event energy spectra calculated using Ajith et al. (2011), as opposed to the frequency domain 𝒫_d quantity calculated from the event waveform.

The reweighting operation is directly implemented in the Monte Carlo evaluation described above, which is valid as long as a sufficient number of samples Θ_j are used:

$$\begin{array}{c}\hfill {\displaystyle \int \mathrm{d}\mathrm{\Theta }}{w}_{0\to 1}(\mathrm{\Theta })p_{\mathrm{d}}(\mathrm{\Theta }|{\mathrm{\Lambda }}_0){\mathcal{P}}_{\mathrm{d}}(\mathrm{\Theta })\approx {\displaystyle \sum _j}{w}_{0\to 1}({\mathrm{\Theta }}_j){\mathcal{P}}_{\mathrm{d}}({\mathrm{\Theta }}_j).\end{array}$$(12)

This allows us to evaluate the (costly) 𝒫_d spectra only once, and reweight the contribution of each wave according to a desired target distribution.

In practice, we rely on the source-frame population distributions to sample the GW parameters. To convert these to detector frame, we evaluate the relevant Jacobian matrix assuming a fixed cosmology,

$$\begin{array}{c}\hfill p_{\mathrm{d}}(\mathrm{\Theta }|\mathrm{\Lambda })=p_{\mathrm{s}}(\mathrm{\Theta }|\mathrm{\Lambda })\frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{s}}}{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{d}}}.\end{array}$$(13)

4.2. Effective sample size and sample variance

As a check of the performance of our reweighting approach, we estimate the effective sample size N_eff for different number of samples N_s and different Ω_GW spectra, and ensure N_eff ≫ 1. The effective sample size is calculated from the weights as N_eff = Σw²/(Σw)². We find, for a fixed reweighting set of Λ hyperparameters, that N_eff ≈ 2 × 10⁴ for N_s = 5 × 10⁴, N_eff ≈ 4 × 10⁴ for N_s = 10⁵, and N_eff ≈ 3 × 10⁵ for N_s = 10⁶. In practice, these numbers depend on the size of the parameter space probed by the reweighting. In all implementations shown in this paper we use N_s = 1 × 10⁶ unless otherwise stated, and we have checked that the order of magnitude of N_eff reported here remains reliable for all results shown.

The Ω_GW spectrum is by definition a stochastic observable, and thus presents an intrinsic sample variance⁷. In particular, ${\mathrm{\Omega }}_{\mathrm{GW}}^{\mathrm{BBH}}$ is dominated by a Poisson process given the low merger rate of black hole binaries. We estimate the intrinsic variance of ${\mathrm{\Omega }}_{\mathrm{GW}}^{\mathrm{BBH}}$ associated with different number of samples (which can be directly converted to different observation times, assuming a total merger rate) by calculating the spectrum using different sample draws from a fixed set of population priors. These are shown in Fig. 1 (left panel) using 10⁵ and 10⁶ samples, where the shading indicates the 95% interval over 1000 sets. In this case, increasing the number of samples by a factor of 10 decreases the width of the 95% interval by 52% on average, for frequencies between 20 and 500 Hz. This specific example corresponds to the following set of hyperparameters for the PLPP mass model: α = 3.5, β = 1, δ_m = 4.5, λ = 0.04, m_max = 100, m_min = 4, m_pp = 34, σ_pp = 4. The redshift model is fixed to the default MD model defined above with a local merger rate of R₀ = 15 Gpc⁻³ yr⁻¹. In Fig. 1 we further compare these sample variance uncertainty bands to the 95% confidence on the spectrum including uncertainty on the local merger rate R₀, drawn from Abbott et al. (2023a). Unsurprisingly, when using these large sample numbers, the sample variance is subdominant compared to population parameter uncertainty.

We compare this intrinsic uncertainty to reweighting. As can be seen in Fig. 1 (right panel), reweighted curves for ${\mathrm{\Omega }}_{\mathrm{GW}}^{\mathrm{BBH}}$ for the given set of Λ hyperparameters lie within the 95% sample uncertainty on the spectrum, for different values of N_s. This implies the reweighted ${\mathrm{\Omega }}_{\mathrm{GW}}^{\mathrm{BBH}}$ spectrum for a given population model is within the intrinsic error on that spectrum, and is thus a fair approximation to make.

As we focus on BBHs in this paper, we drop the BBH subscript from ${\mathrm{\Omega }}_{\mathrm{GW}}^{\mathrm{BBH}}$ in what follows. We refer to the BBH spectrum unless otherwise specified.

5. Background projections using popstock

We studied the dependence of the amplitude, spectral shape, and uncertainty of Ω_GW on various models and data products. These studies will fundamentally inform compact binary population parameter estimation campaigns with upcoming GW datasets, for example in the style of Callister et al. (2020) and Abbott et al. (2021c), which use constraints on (or, in the future, measurements of) the Ω_GW spectrum.

We considered here a frequency range of 10 − 2000 Hz as this corresponds to the sensitivity of second-generation ground-based GW detectors, such as the current configurations of the LIGO, Virgo, and KAGRA instruments as well as their near-future improvements.

5.1. Mass and redshift models

With popstock we were able to rapidly assess the impact of different mass and redshift models on the projected Ω_GW. Here we compared the PLPP and BPL mass models, and the UICV and MD redshift models, all introduced in Sect. 3. When comparing mass models, we fixed the mass model hyperparameters, while including the uncertainty on the local merger rate R₀ from Abbott et al. (2023a) ⁸ assuming the MD redshift evolution with all other parameters fixed to the values discussed above. When comparing redshift models, we instead fixed the redshift model hyperparameters, while including the uncertainty from Abbott et al. (2023a) on the PLPP mass model. We deliberately chose values for certain Λ hyperparameters that are unrealistic and not favoured by the current data to showcase the effect different mass and redshift models may have on the Ω_GW spectral shape.

A comparison between the PLPP and BPL mass models is shown in Fig. 2. The PLPP parameters are fixed to the same set used in Sect. 4.2, and BPL to α₁ = −2, α₂ = −1.4, β = 1, b = 0.4 (see Sect. 3 for details on the parameters). As PLPP is commonly used as a mass model when generating Ω_GW forecasts (as in Abbott et al. 2021c, 2023a), we took this as the fiducial model to produce Ω_GW spectra and compared those obtained with the BPL mass model against these. We note that in particular the choice of α₂ >  α₁ for BPL here implies a larger amount of high-mass binaries in the distribution, as seen in the left panel of Fig. 2, where the primary mass probability distributions are shown. These more massive binaries merge at lower frequencies and their emission is further redshifted into the lower end of the frequency range considered here. We recall that, for example, we detect an equal-mass binary with true component masses of 70 M_⊙ merging around ∼60 Hz at z = 0 and ∼30 Hz at z = 1 (see e.g. Renzini et al. 2024 for more considerations along this line). As seen in the right panel of Fig. 2, this both boosts the amplitude of Ω_GW at all frequencies below a few hundred Hz, and changes the spectral shape at these frequencies, when compared to the PLPP mass spectrum. Specifically, the typical turnover in the spectrum corresponding to the frequency at which most binaries have merged is broken into two turnovers: one for the higher mass binaries (below 100 Hz) and one for the lower mass ones (around 300 Hz). This effect is certainly fuelled by the unrealistic parameter choice made for BPL (α₂ >  α₁). Comparatively, the PLPP model gives rise to a single turnover with a plateau between ∼100 and 400 Hz, which presents small features (‘wiggles’), which are related to the redshifting of the peak at 33 M_⊙. Futhermore, the spectral index at lower frequencies is more peaked than that of the PLPP model. A simple broken-power-law fit to the two sets of curves shown in the right panel of Fig. 2 yields α = 0.61 and α = 0.76 for the lower frequency region of the Ω_GW spectrum corresponding to the PLPP and BPL mass models, respectively. A broader discussion on data-informed spectral indices is postponed to Sect. 5.3.

Fig. 2. Impact of the primary mass distribution on the ΩGW spectrum. The left panel shows the two primary mass model probability densities used throughout. The right panel shows the 95% confidence intervals for ΩGW using the two mass models, including uncertainty on the local merger rate from Abbott et al. (2023a).

The uncertainty on the local rate R₀ implies that there is significant overlap between the 95% credible envelope of the spectrum from these two mass models, suggesting it would be challenging to distinguish mass spectrum features from redshift ones from a measurement of Ω_GW alone. The overlap would be even greater when including the full uncertainty on the redshift model parameters. However, if a large amplitude and large spectral index (i.e. α >  2/3) Ω_GW is observed at low frequencies, we expect a mass model that admits high mass binaries (such as the BPL one showed here) to be favoured.

In Fig. 3 we compare the effect of the UICV and MD redshift models on Ω_GW. We fix the local merger rate to R₀ = 15 Gpc⁻³ yr⁻¹, and compare a UICV rate evolution to the default MD evolution (see Sect. 3), while we include the uncertainty on the PLPP mass model from Abbott et al. (2023a). Most notably, the UICV model impacts the overall amplitude of the Ω_GW spectrum across all frequencies. In this test case, the decrease in amplitude when assuming UICV is approximately constant (and equal to a factor of ∼4) between 10 and 100 Hz, and is due to the lower merger rate between 1 <  z <  4. This effect is much greater than the impact of the uncertainty on the mass model, confirming that a measurement of Ω_GW provides significant information on the merger rate redshift evolution (as also seen in Callister et al. 2016, 2020; Renzini et al. 2024). The turnover in the spectrum appears shifted to slightly higher frequencies in UICV, possibly due to the slightly higher rate fraction at low redshift compared to MD, which instead increases as a power law⁹. Otherwise, the different redshift models do not appear to cause large variations in the overall shape of the spectrum, suggesting that the mass spectrum dominates these features.

Fig. 3. Impact of the merger rate redshift distribution model on the ΩGW spectrum. The left panel shows the two redshift evolution model probability densities used throughout. The right panel shows the 95% confidence intervals for ΩGW using the two merger rate models, including uncertainty on the PLPP mass model from Abbott et al. (2023a), as described in Sect. 3.

A natural extension of this study is to consider BBH mass spectra that evolve with redshift. This will mix the effects seen here when considering independent contributions, and in principle will need to be appropriately included in parameter estimation studies to avoid biases.

5.2. Waveform approximants

The choice of waveform approximant, while central in certain studies of individual compact binary merger events, has been explored very little in the context of the compact binary background signal. In previous work (see e.g. the approximations made in Meacher et al. 2015; Callister et al. 2016; Mukherjee & Silk 2021), it was deemed sufficient to capture solely the evolution of the GW amplitude as a function of frequency as the background Ω_GW contains no phase information. This evolution can be tracked analytically up to arbitrary post-Newtonian (PN) order, considerably speeding up the calculation of Ω_GW compared to calculating full waveforms for large sets of events. Specifically, most works employ the amplitude component of frequency-domain inspiral-merger-ringdown (IMR) waveforms (Ajith et al. 2008, 2011) defined analytically by parts, where the transition of the GW from one phase to the next is set by the specific intrinsic parameters of the binary (mass and spin). As the background has remained a weak signal in the current LVK data, a precise quantification of the systematic differences between background estimates with different waveform approximants has not been necessary. However, as detector sensitivities improve and detection becomes a real possibility, all modelling systematics are important to quantify (see also Zhou et al. 2023; Song et al. 2024). Here, we investigated the effect of the waveform approximant using popstock and confirmed whether it is subdominant to the impact of population uncertainties.

We focus here on the IMRPhenom family of waveforms commonly used in the literature to compute Ω_GW, as well as an effective-one-body (EOB) numerical-relativity (NR) waveform model. In all cases, we omitted spin effects, setting both black hole spins to 0. The following specific waveforms were used:

• IMRPhenomA (Ajith et al. 2008): The first IMR waveform, developed for GW data analysis in the frequency domain for non-spinning binaries. Here the amplitude is expanded to leading post-Newtonian (0-PN) order, implying that the inspiral phase presents the characteristic f^2/3 trend (in Ω_GW units).

• IMRPhenomB (Ajith et al. 2011): A direct successor to IMRPhenomA, this waveform includes higher order corrections in the amplitude term up to 1.5-PN and includes non-zero aligned spin. These correct the spectral shape of the waveform amplitude, as a function of both mass and spin.

• IMRPhenomD (Khan et al. 2016): A recent waveform including corrections up to 3-PN order in the amplitude and a more sophisticated fit to numerical relativity compared to IMRPhenomA and B.

• SEOBNRv2 (Pürrer 2016): An EOB NR waveform for spin-aligned BBHs, calculated numerically in the frequency domain.

A comparison between the Ω_GW spectra calculated for the same BBH population using the waveforms above is presented in Fig. 4. We show the 95% confidence contours on the population Ω_GW including the uncertainty on the mass model (PLPP) and the local merger rate R₀ assuming a fixed MD redshift evolution (as defined above).

Fig. 4. Impact of the waveform model on the ΩGW spectrum. The shading indicates the 95% confidence on the spectrum including the uncertainty on the PLPP mass model and the local merger rate, assuming a fixed MD redshift evolution.

We find that when employing the 0-PN IMRPhenomA waveform, the background signal is overestimated at all frequencies by up to 50% in the range 10 <  f <  1000 Hz compared to IMRPhenomB; this is due to the missing amplitude corrections to the inspiral phase (see e.g. the differences in Eq. (4.13) of Ajith et al. 2008 and Eq. (1) of Ajith et al. 2011). While the amplitude estimates for IMRPhenomA and B agree at f ≡ 0, they diverge for f >  0 as the amplitude evolution with frequency is Ω_A ∝ f^2/3 for IMRPhenomA and Ω_B ∝ f^{α <  2/3} for IMRPhenomB. The trend of α depends on the specific mass, redshift, and spin realisation, as discussed in Sect. 5.1. This result shows that the somewhat basic assumption that the compact binary background at frequencies under ∼100 Hz is well-approximated by a two-thirds power law can be upgraded and informed by likely population models to optimise background searches.

Differences due to the use of IMRphenomB/D and SEOBNRv2 approximants are comparable to each other and would be hard to distinguish from population uncertainties. Nevertheless, we comment that the different NR calibration used in IMRphenomD compared to B is evident in the impact due to the inspiral phase on the GWB signal, as the frequency evolution at low frequency is slightly modified, and SEOBNRv2 estimates an overall lower background than the IMRPhenom waveforms.

The impact of including higher-order modes in the waveform calculation on the background spectrum was found to be negligible; a comparison between spectra calculated with the IMRPhenomD and IMRPhenomXPHM waveforms is included for completeness in Appendix B.

5.3. O3 population samples

We conclude our analyses by combining the uncertainties on the mass and redshift distributions drawn directly from the LVK GWTC-3 population analysis (Abbott et al. 2023a). We limit our focus to the MD redshift model for BBHs as this is the only model we have viable posterior samples for. In the case of the low-redshift merger rate parameters (R₀, γ), we used the samples from the power-law redshift inference results released in LIGO Scientific Collaboration, Virgo Collaboration, KAGRA Collaboration (2021) for the power-law redshift model, while when including high-redshift features (z_max, κ) we used the results obtained performing inference on the entire MD model, as done in Abbott et al. (2021c).

The results obtained progressively varying the redshift hyperparameters are shown in Fig. 5 (left panel). We included uncertainty over the entire PLPP mass hyperparameter space, while varying the following:

Fig. 5. Uncertainty on the expected ΩGW spectrum from BBHs due to uncertainty on the merger rate evolution parameters. Left: 95% confidence levels on the projected ΩGW spectrum including uncertainty on the PLPP mass model and assuming a MD merger rate model with different levels of uncertainty. The hatched outline reports previous results published in Abbott et al. (2023a). Right: zoomed-in comparison at low frequencies of the uncertainty on the ΩGW spectrum when varying R0 vs. R0 and γ, reporting average spectral indices referred to these two contours.

• (i)

  just the local merger rate R₀, assuming a fixed broken-power-law merger rate evolution with parameters fixed to those discussed in Sect. 3;

• (ii)

  both R₀ and the local power-law index γ, keeping the high-redshift parameters fixed;

• (iii)

  all parameters for the merger rate, including the turnover redshift z_peak and high-redshift power-law index κ.

In cases (i) and (ii) the samples were drawn from the power-law redshift model posteriors of Abbott et al. (2023a). In case (iii), the samples were instead drawn from a full GWTC-3, O1–O2–O3 joint stochastic-population analysis (similar to Abbott et al. 2021c; for details on how this analysis was carried out, see Callister et al. 2020) as these also include posteriors on the higher redshift evolution of R(z) (Callister 2024).

We compared the 95% confidence levels on Ω_GW in case (i) with published results (shown in Fig. 5 in hatched outline¹⁰) and find them to be consistent. We note that the corresponding LVK contours draw from the same PLPP mass posterior and local merger rate posterior, but assume a different redshift evolution (see the original paper discussion Abbott et al. 2023a), which explains the small differences between the curves at low frequency and the different turnover trend at high frequency, which is dominated by low-redshift effects. Specifically, the redshift model used sampled a time-delay distribution between binary formation and merger, where the binary formation rate was fixed to follow the star formation rate of Vangioni et al. (2015), and the time-delay distribution was in the shape of $p(t_{\mathrm{d}})\propto t_{\mathrm{d}}^{-1}$. This model cannot be well approximated with a simple broken power law. The LVK contour was calculated assuming the IMRPhenomB analytic waveform model. Case (ii) and case (iii) give almost identical contours, which implies that the population analysis Abbott et al. (2023a) and the stochastic constraints Abbott et al. (2021c) carry little information about the high-redshift evolution of the merger rate. Furthermore, they show that the uncertainty on the local merger rate evolution alone could account for an increase of up to a factor of ∼5 in the Ω_GW spectrum. This could have considerable consequences on the detectability of the signal. In the right panel of Fig. 5 we show a zoomed-in image of the [20, 200] Hz portion of the spectra in cases (i) and (ii), and provide results of single power-law fits to the envelope of Ω_GW curves. The average α spectral indices found are consistent with each other: α_(i) = 0.59 ± 0.02 for (i), and α_(ii) = 0.60 ± 0.03 for (ii). This confirms that γ has no impact on the low-frequency spectral shape of Ω_GW, but only on its amplitude.

We repeated the Ω_GW calculation varying R₀ and γ using the BPL mass model, sampling over the joint mass and redshift posteriors obtained by performing population inference on the GWTC-3 catalogue (samples are publicly available as a part of the example sample sets in the popstock package repository). We compared these results with the PLPP results described above in Fig. 6, overlaying the 2σ LVK power-law-integrated sensitivity curves (Thrane & Romano 2013) already shown in Abbott et al. (2023a). These track the present and future upgrades to the LIGO and Virgo facilities, where ‘O3 sensitivity’ is given by the O3 measured spectra, ‘Design HLV’ is produced assuming projections shown in Abbott et al. (2018) and is expected to approximate the sensitivity at the end of the O4 observing run (currently ongoing), and ‘Design A+’ refers to the sensitivity projected for the next observing run O5 assuming one year of continuous data and all planned improvements to the detectors are implemented and successful.

Fig. 6. Projections of the background ΩGW spectrum, given our knowledge of the compact binary population. The shaded regions outline 90% credible bands for the GWB from BBHs and BNS (in pink), including the uncertainty on the mass and redshift models for these sources using the samples released in Abbott et al. (2023a), LIGO Scientific Collaboration, Virgo Collaboration, KAGRA Collaboration (2021). For BBHs we report the uncertainty due to two mass models: the PLPP mass model, assuming a MD redshift model with uncertain local merger rate R0 (orange), and also uncertain low-redshift power-law index γ (green); and the BPL mass model, assuming a MD redshift model with uncertain R0 and γ. Current and projected sensitivity curves are included for reference.

The BPL mass model predicts a systematically larger background, by a factor of 1.4 on average, which hints at the possibility of a louder signal than previously projected, and thus the prospect of a detection of the stochastic signal before reaching the Design A+ LIGO-Virgo sensitivity. For completeness, we also include the expected background signal arising from binary neutron stars (BNS) in Fig. 6. This signal is strongly dominated by uncertainties given the very few detections of BNS mergers (Abbott et al. 2017, 2020). The projection is calculated employing the NRTidalv2 model discussed in Dietrich et al. (2019). We assume a uniform mass model between 1 and 2.5 M_⊙, as in Abbott et al. (2023a), and a merger rate model corresponding to a time-delayed SFR as used for projections presented in Abbott et al. (2021c) ¹¹. We draw the local merger rate R₀ from the corresponding posterior samples presented in Abbott et al. (2023a), LIGO Scientific Collaboration, Virgo Collaboration, KAGRA Collaboration (2021); for reference, we refer to the samples that set the $R_0=105.5_{-83.9}^{+190.2}$ Gpc⁻³ yr⁻¹ upper limit.

6. Conclusions

We present a novel method and code-base for rapid calculations of the background spectrum of inspiralling and coalescing compact binaries, starting from a given population model and hyperparameter sets.

We quantified the joint uncertainty on the Ω_GW spectrum from both the mass and redshift distributions of the BBH population, given the most recent results from the LVK collaboration. Predictably, the uncertainty on the local merger rate and its evolution (together with the uncertainty on the mass model) dominated the expected amplitude of the spectrum, and can have significant implications on detectability.

Furthermore, we find that for the preferred mass model (power-law-plus-peak), the low-frequency spectral index of the stochastic background signal is α = 0.6. Previous detection approaches assumed α = 2/3 (for example Abbott et al. 2021b,c; Renzini & Contaldi 2019); this result was based on the waveform used to calculate the expected GWB and its PN order expansion, as first seen in Phinney (2001) and then Sesana et al. (2008). We find that, when employing 0-PN order waveforms, there is a tension between projections for Ω_GW from the presently observed population which competes with population uncertainty itself. The mismatch between the treatment of the late inspiral phase in IMR waveforms is particularly significant as it is present throughout the entire Ω_GW spectrum and in particular at lower frequencies, where current detector sensitivities peak.

Differences between the specific background realisations and the number of samples also produce different projections that may give rise to small tensions in the higher frequency range, where the spectra present a turnover that is highly dependent on the binary mass distribution and local features of the merger rate. Current-generation detectors are not sensitive to this region of the spectrum; however, this will have significant implications for next-generation interferometers.

In conclusion, the specific population properties of the coalescing compact binary population as well as the specific realisation during our observations will play a role in detection capabilities. In particular, within current binary black hole population uncertainties, a low-redshift amplification of the merger rate and a larger population of higher-mass binaries contribute to a significant boost in the background amplitude, in the LVK sensitivity band, which could lead to early detection. On the other hand, more astrophysically motivated BBH rate evolution models relate the merger rate to binary progenitor features, and reparametrise the merger rate density in terms of, for example, the time-to-merger delay distribution and the host galaxy metallicity (Fishbach & Kalogera 2021; Chruślińska 2024). These models have recently been employed in joint analyses of the GWTC-3 catalogue and LVK stochastic upper limits (Turbang et al. 2024), and may provide alternative forecasts of the uncertainty on the Ω_GW spectrum, as we gather more GW data. These models will progressively be included and updated in popstock. With popstock, we provide the GW community engaged in GW source modelling, data analysis, and astrophysical interpretation with a user-friendly tool for rapid background spectrum evaluation and easy integration of new models as our understanding of the GW universe expands.

Acknowledgments

We thank Patrick Meyers and Alan Weinstein for invaluable discussions and insight. We thank Thomas Callister for providing the population samples in Callister (2024), and for carefully reading our work. We thank Nicholas Loutrel for consulting on waveform models and their features. AIR is supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 101064542, and acknowledges support from the NSF award PHY-1912594. JG is supported by NSF award No. 2207758. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation.

References

Appendix A: Deriving the energy in GWs

We derive here the energy spectrum dE/df carried by gravitational waves, in vacuum (see e.g. Maggiore 2007; Saffer et al. 2018; Carroll 2004). We start by expanding the perturbed metric g_μν to second order,

$$\begin{array}{c}\hfill g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }^{(1)}+h_{\mu \nu }^{(2)},\end{array}$$(A.1)

where η_μν is the Minkowski flat metric and it is assumed that the perturbation $h_{\mu \nu }^{(i)}$ is ith order in some small parameter controlling the scale of h_μν. Substituting the above into the Einstein field equation gives

$$\begin{array}{c}\hfill G_{\mu \nu }\left[h^{(1)}\right]+G_{\mu \nu }\left[{\left(h^{(1)}\right)}^2\right]+G_{\mu \nu }\left[h^{(2)}\right]=8\pi G{\tau }_{\mu \nu },\end{array}$$(A.2)

Where the first term is the Einstein tensor linear in the first order perturbation, the second is the Einstein tensor terms quadratic in the first order perturbation, and the third term is linear in the second-order perturbation. In vacuum, τ_μν = 0, and the solution for h⁽¹⁾ (i.e. the plane wave solution) reduces the first term to $G_{\mu \nu }\left[h^{(1)}\right]=\square h_{\mu \nu }^{(1)}=0$, in the Lorenz gauge. We can therefore rearrange the above into a form that resembles the Einstein field equation,

$$\begin{array}{c}\hfill G_{\mu \nu }\left[h^{(2)}\right]=-G_{\mu \nu }\left[{\left(h^{(1)}\right)}^2\right],\end{array}$$(A.3)

where the first-order term $h_{\mu \nu }^{(1)}$ squared effectively forms a stress-energy (RHS) that sources the second order curvature (LHS). In this analogy to the RHS of the Einstein equation, we can define the effective stress energy (pseudo)-tensor of GWs Isaacson (1968):

$$\begin{array}{c}\hfill {\tau }_{\mu \nu }\equiv -\frac{1}{8\pi G}{G}_{\mu \nu }\left[{\left(h^{(1)}\right)}^2\right].\end{array}$$(A.4)

A nice feature here is that the LHS of Eq. A.3 satisfies the contracted Bianchi identities and therefore the RHS is divergence-free and can be interpreted as conserving energy according to some observer. Expanding out (A.4) gives

$$\begin{array}{c}\hfill {\tau }_{\mu \nu }=\frac{c^4}{32\pi G}\langle {\partial }_{\mu }{h}_{\alpha \beta }{\partial }_{\nu }{h}^{\alpha \beta }\rangle .\end{array}$$(A.5)

The conservation law ∂_μt^μν = 0 implies

$$\begin{array}{c}\hfill {\partial }_0{\tau }^{00}+{\partial }_i{\tau }^{i0}=0,\end{array}$$(A.6)

where τ⁰⁰ can be interpreted as a volumetric energy density, such that the energy E is defined as E = ∫d³xτ⁰⁰ and therefore the associated power is

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dt}}={\partial }_0{\displaystyle \int d^3}x{\tau }^{00}.\end{array}$$(A.7)

Substituting into Eq. A.6 yields

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dt}}+{\displaystyle \int d^3}x{\partial }_i{\tau }^{i0}=0,\end{array}$$(A.8)

which simplifies to

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dt}}+{\displaystyle \int d^3}x{\partial }_z{\tau }^{z0}=0\end{array}$$(A.9)

for a wave moving along direction $\widehat{z}$. We employ the divergence theorem to convert the volume integral into a surface integral,

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dt}}+\widehat{z}{\displaystyle \int dA}{\tau }^{z0}=0,\end{array}$$(A.10)

which gives an expression for the energy flux (energy per unit time per unit area):

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dAdt}}=-{\tau }^{z0}\widehat{z}\equiv -{\tau }^{00}\widehat{z},\end{array}$$(A.11)

as ∂₀h_ij = −∂_zh_ij = −∂⁰h_ij holds for a wave solution. Considering our gauge, the only surviving terms are those for μ, ν = 1 or 2,

$$\begin{array}{cc}\hfill {\tau }^{00}=\frac{c^4}{32\pi G}\langle {\partial }_0{h}_{11}{\partial }_0{h}^{11}+{\partial }_0{h}_{12}{\partial }_0{h}^{12}+{\partial }_0& h_{21}{\partial }_0{h}^{21}+\hfill \\ \hfill & +{\partial }_0{h}_{22}{\partial }_0{h}^{22}\rangle ,\hfill \end{array}$$(A.12)

and substituting in the polarisation components of the wave h_μν in the TT gauge yields

$$\begin{array}{c}\hfill {\tau }^{00}=\frac{c^2}{16\pi G}\langle |{\dot{h}}_{+}{|}^2+{|{\dot{h}}_{\times }|}^2\rangle .\end{array}$$(A.13)

Solving for the energy flux of gravitational waves of Eq. (A.10) gives

$$\begin{array}{c}\hfill \left|\frac{\mathit{dE}}{\mathit{dAdt}}\right|=-\frac{c^3}{16\pi G}\langle |{\dot{h}}_{+}{|}^2+{|{\dot{h}}_{\times }|}^2\rangle ,\end{array}$$(A.14)

where we have switched to the absolute value of this quantity with the understanding that GWs are removing energy from the system. The surface area energy density is then defined as

$$\begin{array}{c}\hfill \frac{\mathit{dE}}{\mathit{dA}}={\displaystyle \int dt}\frac{c^3}{16\pi G}\langle {\dot{h}}_{+}{(t)}^2+{\dot{h}}_{\times }{(t)}^2\rangle ,\end{array}$$(A.15)

To expand the above we recall that, for a plane wave,

$$\begin{array}{c}\hfill \dot{\stackrel{\sim }{h}}(f)={\displaystyle \int dt}\dot{h}(t)e^{-i\omega t}=i\omega {\displaystyle \int dt}h(t)e^{-i\omega t}=i\omega \stackrel{\sim }{h}(f),\end{array}$$(A.16)

where in the final equivalence we have directly employed the definitaion of a Fourier transform. Recalling Parseval’s theorem, we can write the surface area energy density in terms of the Fourier transform $\stackrel{\sim }{h}(f)$,

$$\begin{array}{cc}\hfill \frac{\mathit{dE}}{\mathit{dA}}& ={\displaystyle {\int }_{-\infty }^{\infty }}df\frac{c^3{\omega }^2}{16\pi G}\langle {\stackrel{\sim }{h}}_{+}{(f)}^2+{\stackrel{\sim }{h}}_{\times }{(f)}^2\rangle \hfill \\ \hfill & =\frac{\pi c^3}{2G}{\displaystyle {\int }_0^{\infty }}df{f}^2\langle {\stackrel{\sim }{h}}_{+}{(f)}^2+{\stackrel{\sim }{h}}_{\times }{(f)}^2\rangle ,\hfill \end{array}$$(A.17)

where the integral is taken over the sphere surrounding the source. Note that h₊ and h_× terms include dependence on inclination ι and reference phase ϕ₀ (or, equivalently, the observer’s position along the azimuth) and therefore must be included in the integral over the area.

Appendix B: Impact of higher order modes on the GWB spectrum

For the sake of completeness, we append here findings on the impact on the Ω_GW spectrum due to the inclusion of higher order modes in the waveform model. Higher order modes are subdominant harmonics excited during GW emission, where the dominant harmonic is the ℓ = 2, m = 2 mode London et al. (2018). We compare the 95% confidence bands shown in 5.2 for the IMRPhenomD waveform model with bands obtained using the IMRPhenomXPHM approximant Pratten et al. (2021), which includes the (ℓ, |m|) = (2,  2), (2,  1), (3,  3), (3,  2), (4,  4) modes. As may be seen in Fig. A.1, for the population of binaries considered in 5.2 which in particular is non-spinning and non-precessing, there are negligible differences between the use of IMRPhenomD and IMRPhenomXPHM. This is particularly evident when comparing the Ω_GW spectrum calculated using the same event samples employing the two waveforms, shown as dashed and dotted curves in the left panel of Fig. A.1. The percent difference between these two curves is reported in the right panel of Fig. A.1, which remains consistently below 10% across the spectrum and is under 3% for frequencies below 100 Hz.

Fig. A.1. Impact of the inclusion of higher order modes in the waveform model employed to evaluate the ΩGW spectrum. On the left: 95% confidence on the spectrum including the uncertainty on the PLPP mass model and the local merger rate, assuming a fixed MD redshift evolution. On the right: percent difference %ΔΩGW(f) between ΩGW spectra calculated using the same event samples, shown as dashed and dotted curves on the left panel.

All Tables

All Figures

Fig. 1. Impact of sample variance and reweighting on the ΩGW spectrum. Left: 95% confidence on the spectrum calculated using 105 (106) samples in yellow (green) drawn from a fixed hyperparameter distribution, compared to the 95% confidence on the spectrum including uncertainty on the local merger rate parameter in brown. Unsurprisingly, the sample variance varies significantly with number of samples and when using large sample numbers is subdominant compared to population parameter uncertainty. Right: sample variance from the left panel compared to reweighted estimates of the ΩGW spectrum. The reweighted spectra lie neatly within the sample variance uncertainty bounds, implying that the reweighted spectrum is indistinguishable from a regularly sampled spectrum, with these sample numbers.

In the text

	Fig. 2. Impact of the primary mass distribution on the ΩGW spectrum. The left panel shows the two primary mass model probability densities used throughout. The right panel shows the 95% confidence intervals for ΩGW using the two mass models, including uncertainty on the local merger rate from Abbott et al. (2023a).
In the text

	Fig. 3. Impact of the merger rate redshift distribution model on the ΩGW spectrum. The left panel shows the two redshift evolution model probability densities used throughout. The right panel shows the 95% confidence intervals for ΩGW using the two merger rate models, including uncertainty on the PLPP mass model from Abbott et al. (2023a), as described in Sect. 3.
In the text

	Fig. 4. Impact of the waveform model on the ΩGW spectrum. The shading indicates the 95% confidence on the spectrum including the uncertainty on the PLPP mass model and the local merger rate, assuming a fixed MD redshift evolution.
In the text

Fig. 5. Uncertainty on the expected ΩGW spectrum from BBHs due to uncertainty on the merger rate evolution parameters. Left: 95% confidence levels on the projected ΩGW spectrum including uncertainty on the PLPP mass model and assuming a MD merger rate model with different levels of uncertainty. The hatched outline reports previous results published in Abbott et al. (2023a). Right: zoomed-in comparison at low frequencies of the uncertainty on the ΩGW spectrum when varying R0 vs. R0 and γ, reporting average spectral indices referred to these two contours.

In the text

Fig. 6. Projections of the background ΩGW spectrum, given our knowledge of the compact binary population. The shaded regions outline 90% credible bands for the GWB from BBHs and BNS (in pink), including the uncertainty on the mass and redshift models for these sources using the samples released in Abbott et al. (2023a), LIGO Scientific Collaboration, Virgo Collaboration, KAGRA Collaboration (2021). For BBHs we report the uncertainty due to two mass models: the PLPP mass model, assuming a MD redshift model with uncertain local merger rate R0 (orange), and also uncertain low-redshift power-law index γ (green); and the BPL mass model, assuming a MD redshift model with uncertain R0 and γ. Current and projected sensitivity curves are included for reference.

In the text

Fig. A.1. Impact of the inclusion of higher order modes in the waveform model employed to evaluate the ΩGW spectrum. On the left: 95% confidence on the spectrum including the uncertainty on the PLPP mass model and the local merger rate, assuming a fixed MD redshift evolution. On the right: percent difference %ΔΩGW(f) between ΩGW spectra calculated using the same event samples, shown as dashed and dotted curves on the left panel.

In the text