© The Authors 2025

1. Introduction

Gravitational waves (GWs) from compact binary coalescences (CBCs) offer a unique opportunity to study the expansion of the Universe, especially in light of ongoing tensions in cosmology (Di Valentino et al. 2021). Unlike electromagnetic (EM) observables, CBCs detected through GWs carry direct information about their luminosity distance; as such, they are referred to as standard sirens (Schutz 1986; Holz & Hughes 2005). When combined with redshift information, GW data can therefore be used to probe the luminosity distance–redshift relation – and hence key cosmological parameters such as the Hubble constant, H₀. However, GWs alone do not provide the source redshift. Various methods have been developed to estimate redshifts for GW events, which are essential for this purpose.

Bright sirens are a special class of standard sirens in which an EM counterpart is observed along with the GW emission. This was the case of the binary neutron star merger GW170817 detection in association with a gamma-ray burst and a subsequent kilonova (Abbott et al. 2017a). The host galaxy of this source was identified, allowing us to obtain a precise redshift measurement and an H₀ estimate of ${70}_{-8}^{+19}$ km/s/Mpc (Abbott et al. 2017b). So far, GW170817 is the only known event with an associated EM counterpart.

Most often, other detections involve dark sirens, such as binary black hole (BBH) mergers, with no EM counterparts. For these types of events, several methods have been proposed to obtain a statistical redshift estimation for the source. One approach involves cross-matching the event’s sky localization with galaxy catalogs, assigning probabilities based on the distribution of potential host galaxies (Abbott et al. 2017c, 2021a, 2023a; Soares-Santos et al. 2019; Gray et al. 2020, 2022, 2023; Palmese et al. 2020; Finke et al. 2021; Mukherjee et al. 2024; Mastrogiovanni et al. 2023; Palmese et al. 2023; Bom et al. 2024). Clustering-based techniques, another method, exploit correlations between GW events and large-scale structure to constrain redshifts (Namikawa et al. 2016a,b; Oguri 2016; Zhang 2018; Scelfo et al. 2020; Bera et al. 2020; Libanore et al. 2021; Mukherjee et al. 2021; Cigarrán Díaz & Mukherjee 2022; Ferri et al. 2025; Ghosh et al. 2025; Zazzera et al. 2025).

In this paper we focus on the spectral sirens method (Taylor et al. 2012; Farr et al. 2019; Mastrogiovanni et al. 2021; Mancarella et al. 2022, 2024; Ezquiaga & Holz 2021, 2022; Mali & Essick 2025). This approach uses the inherent degeneracy between a source’s redshift, z, and its detector frame mass, m_det, as well as its source frame mass, m_s, which is linked by the relation m_det = (1 + z)m_s. The sources’ redshifts are indirectly inferred by measuring the detector frame masses and simultaneously modeling the source mass distribution and the cosmological parameters. The effectiveness of this method depends on the choice of phenomenological models used to represent the CBC source mass distribution and their ability to accurately capture the true characteristics of the CBC population (Mukherjee 2022; Gray et al. 2022; Karathanasis et al. 2023; Pierra et al. 2024; Mastrogiovanni et al. 2023; Agarwal et al. 2025). Recent developments also include nonparametric (Farah et al. 2025; Magaña Hernandez & Ray 2024) and multi-subpopulation (Li et al. 2024; Tong et al. 2025) approaches to mass function inference, offering greater flexibility and robustness in capturing the true mass distribution. This approach was applied to measure H₀ using 47 GW sources from the Third LIGO and Virgo Gravitational-Wave Transient Catalog (GWTC–3); when combined with the H₀ measurement from GW170817 and its EM counterpart, $H_0={68}_{-7}^{+12}$ km/s/Mpc (68% credible interval) was obtained (Abbott et al. 2023b).

In addition to resolved sources, we expect to detect a stochastic gravitational-wave background (SGWB) from unresolved sources. The SGWB depends on the properties of the population and the underlying population of CBCs (Romano & Cornish 2017; Christensen 1992; van Remortel et al. 2023). The idea is to explore whether the detection or non-detection of this stochastic background, combined with the information from resolved sources, can help constrain the population properties more accurately and, in turn, refine the measurements of the cosmological parameters (Abbott et al. 2021b). Moreover, as recognized in Cousins et al. (2025), the SGWB also has a direct dependence on H₀ due to the fact that the GW energy density is related to the comoving volume. Studying the interplay between GW signals from CBCs and the SGWB presents an exciting opportunity to improve our understanding of fundamental cosmological and astrophysical parameters.

The method of synthesizing data from individual CBC events with upper limits on the SGWB was initially examined by Callister et al. (2020) and Turbang et al. (2024) in the context of astrophysical properties for CBCs and applied in Ezquiaga (2021) and Zhong et al. (2025). More recently, it was used by Cousins et al. (2025) for a joint cosmological and population analysis with GWTC-3. In our study, we focus more on understanding what population properties and cosmological parameters could be constrained by the SGWB for future observing runs.

This paper is organized as follows. Section 2 outlines the methodology for our integrated analysis of individual BBH mergers and the stochastic background. Section 3 presents how we simulated the data. We discuss our projection studies in Sect. 4 and a reanalysis with this methodology of GWTC-3 in Sect. 5. Finally, in Sect. 6 we summarize our conclusions.

2. Analysis method

We aimed to infer the cosmological and population parameters including resolved and unresolved CBC. We built on the framework established by Callister et al. (2020) and Turbang et al. (2024). We describe the GW dataset as a set, {x}, of N_obs resolved GW detections found in data chunks {x}={x₁, …, x_Nobs} and a set of elements $\widehat{C}(f)$ describing the amount of correlated noise across GW detectors that are obtained for an observing time T_obs. Let Φ denote the hyperparameters that characterize the CBC population and cosmology (e.g., H₀ and Ω_m). We proceeded under the assumption that the joint likelihood of these observations can be factorized as

$$\begin{array}{c}\hfill \mathcal{L}(\{x\},\widehat{C}|\mathrm{\Phi })={\mathcal{L}}_{\mathrm{CBC}}(\{x\}|\mathrm{\Phi }){\mathcal{L}}_{\mathrm{SGWB}}(\widehat{C}|\mathrm{\Phi }),\end{array}$$(1)

where ℒ_CBC represents the likelihood associated with the detection of resolved sources, which, as discussed in Sect. 2.1, follows an inhomogeneous Poisson distribution. The likelihood ℒ_SGWB denotes stochastic likelihood modeled as a Gaussian process, as detailed in Sect. 2.2. We note that the factorization of the likelihood in Eq. (6) assumes that the “stochastic” data are independent of the resolved sources. In other words, we assumed that the amount of stochastic GW signal is not modified after all the resolved sources in the observing period, T_obs, are removed. This might not be the case when the sensitivity of the GW detectors is enough to solve a considerable amount of the population. To address this issue while still using the same likelihood model, we split the GW data collected in an observing time, T_obs, into two disjoint sets with duration T_obs/2 and used one for detecting solved sources and the other for the SGWB.

2.1. Hierarchical Bayesian likelihood

The detection of GWs events is modeled as an inhomogeneous Poisson process with selection biases (Mandel et al. 2019; Vitale et al. 2021). For N_obs detected GWs signals over an observation time, T_obs, the probability of obtaining a specific GWs dataset, {x}, given population hyperparameters Φ is

$$\begin{array}{cc}\hfill {\mathcal{L}}_{\mathrm{CBC}}(\{x\}|\mathrm{\Phi })\propto & e^{-N_{\exp }(\mathrm{\Phi })}{\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}}{T}_{\mathrm{obs}}{\displaystyle \int \mathrm{d}}\theta \mathrm{d}z\hfill \\ \hfill & \times {\mathcal{L}}_{\mathrm{GW}}(x_i|\theta ,z,\mathrm{\Phi })\frac{1}{1+z}\frac{\mathrm{d}{N}_{\mathrm{CBC}}}{\mathrm{d}\theta \mathrm{d}z\mathrm{d}{t}_{\mathrm{s}}}(\mathrm{\Phi }).\hfill \end{array}$$(2)

Here N_exp represents the expected number of CBC detections within the observation time, T_obs. Typically, N_exp is evaluated using a set of Monte Carlo injections in real data and is used to correct for selection biases; more details are provided in Appendix A. The variables z and θ denote the redshift and a set of binary parameters for each CBC, such as the source masses. The individual GW likelihood ℒ_GW(x_i|θ, z, Φ) incorporates uncertainties in the parameters, θ; the final term of the equation corresponds to the CBC rate in the source frame time, t_s. The model that we adopted for the CBC rate is described in more detail in Sect. 3.

2.2. Stochastic gravitational-wave background likelihood

The SGWB is conventionally defined as (Romano & Cornish 2017)

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

where ρ_c = 3H₀²c²/(8πG) is the critical energy density of the Universe, ρ_GW(f) is the SGWB energy density observed at the frequency f, G is Newton’s gravitational constant, c is the speed of light, and H₀ is the Hubble constant. The SGWB energy density modeled from the population of CBC is given by Phinney (2001)

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{GW}}(f)=\frac{f}{{\rho }_{\mathrm{c}}}{\displaystyle {\int }_0^{+\infty }}\mathrm{d}z\frac{R(z)}{H(z)(1+z)}\langle \frac{\mathrm{d}{E}_{\mathrm{s}}}{\mathrm{d}{f}_{\mathrm{s}}}{|}_{f(1+z)}\rangle .\end{array}$$(4)

Here R(z) is the event rate per unit comoving volume and per unit time at redshift, z. The term H(z) is the Hubble parameter that under the assumption of Λ cold dark matter is defined as $H(z)=H_0\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}$ (Wei & Wu 2017). As recognized in Cousins et al. (2025), the SGWB scales as H₀⁻³, this scaling arises from comoving volume density term (∝H₀⁻³). The final term describes the energy spectrum radiated by the CBC population, averaged over the source population at a given redshift. More details on how we calculated the SGWB are provided in Appendix A.

Searches for the SGWB, such as those conducted by the LIGO-Virgo-Kagra (LVK) collaboration (Abbott et al. 2021b), aim to measure the dimensionless energy density, Ω_GW(f), in Eq. (4). Let us label the GW detectors in the network by the index I, the optimal cross-correlation estimator for a baseline IJ is given by (Allen & Romano 1999; Abbott et al. 2021b)

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{IJ}}(f)=\frac{2}{T_{\mathrm{obs}}}\frac{10{\pi }^2}{3H_0^2}\frac{f^3}{{\gamma }_{\mathit{IJ}}(f)}\mathrm{Re}[{\stackrel{\sim }{s}}_I(f){\stackrel{\sim }{s}}_J^{\ast }(f)],\end{array}$$(5)

where T_obs is the observation time, and ${\stackrel{\sim }{s}}_I(f)$ and ${\stackrel{\sim }{s}}_J(f)$ are the Fourier transforms of the data for detectors I and J. The overlap reduction function γ_IJ(f) encodes the baseline geometry of the detector pair IJ (Christensen 1992; Flanagan 1993). The combined cross-correlation estimator using information from all baselines is obtained by

$$\begin{array}{c}\hfill \widehat{C}(f)=\frac{{\displaystyle {\sum }_{\mathit{IJ}}}{\widehat{C}}_{\mathit{IJ}}(f){\sigma }_{\mathit{IJ}}^{-2}}{{\displaystyle {\sum }_{\mathit{IJ}}}{\sigma }_{\mathit{IJ}}^{-2}},\end{array}$$(6)

with an expected value of $⟨\widehat{C}(f)⟩={\mathrm{\Omega }}_{\mathrm{GW}}(f)$. We adopted the shorthand notation ∑_IJ to mean a sum over all independent baselines, IJ. The variance of the estimator for a single baseline, IJ, is

$$\begin{array}{c}\hfill {\sigma }_{\mathit{IJ}}^2(f)\approx \frac{1}{2T_{\mathrm{obs}}\mathrm{\Delta }f}{\left(\frac{10{\pi }^2}{3H_0^2}\right)}^2\frac{f^6}{{\gamma }_{\mathit{IJ}}^2(f)}{P}_I(f)P_J(f),\end{array}$$(7)

under the assumption of a small S/N for the SGWB, where P_I(f) and P_J(f) are the one-sided power spectral densities of detectors I and J, and Δf is the frequency resolution. The combined variance for all pairs, assuming statistical independence, is given by

$$\begin{array}{c}\hfill {\sigma }^2(f)={\left({\displaystyle \sum _{\mathit{IJ}}}\frac{1}{{\sigma }_{\mathit{IJ}}^2(f)}\right)}^{-1}.\end{array}$$(8)

The variance σ²(f) can be used to set the upper limit on the SGWB.

The cross-correlation estimator in Eq. (6) is then used in the likelihood ℒ_SGWB, which is well approximated by a Gaussian distribution (Abbott et al. 2021b):

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{SGWB}}(\widehat{C}|\mathrm{\Phi })\propto exp[-\frac{1}{2}{\displaystyle \sum _k}{\left(\frac{\widehat{C}(f_k)-{\mathrm{\Omega }}_{\mathrm{GW}}(f_k,\mathrm{\Phi })}{\sigma (f_k)}\right)}^2],\end{array}$$(9)

where the sum runs over discrete frequency bins, f_k, and where Ω_GW is our model energy-density spectrum defined in Eq. (4). We note that the stochastic likelihood is not invariant for H₀, but it scales as H₀⁻¹. This scaling is due to the combination of the CBC energy density with the cross-correlation statistic and variance, which also matches the 1/H₀ scaling of the SGWB S/N for detection (Cousins et al. 2025).

3. Framework to simulate GW data for O5 sensitivities

We built two observing scenarios for the simulation that are based on two years of GW data. The first one is called “resolved CBCs” and it consists of 2 years of GW data that are used only to detect resolved sources. The second case is called “combined SGWB and CBCs” and it is generated with 1 year of GW data used to collect resolved sources and 1 year used to detect the SGWB. For these datasets we considered two study cases: a first where we assumed a five-detector network composed of LIGO Livingston, Hanford, India, Virgo, and KAGRA, and a second where we considered a four-detector network composed of LIGO Livingston, Hanford, Virgo, and KAGRA. In the following, we use the five-detector network as our default choice, but we also discuss the possibility of having a four-detector network.

The choice of dividing the GW data into two independent sets for the resolvable sources and the SGWB is driven by the fact that we want to keep the hierarchical likelihood separable. We also have a third validation case with 1 year of GW data employed only for stochastic searches that is used to understand the actual contribution of the SGWB on the population parameters.

The common ingredient needed to simulate data for the solved GW sources and the SGWB is a population and cosmological model for the CBCs. We assumed a Λ cold dark matter cosmology (Wei & Wu 2017) and parametrized the rate of CBCs introduced in Sect. 2 as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{CBC}}}{\mathrm{d}\theta \mathrm{d}z\mathrm{d}{t}_s}=R_0\psi (z,\mathrm{\Phi })p_{\mathrm{pop}}({\mathit{m}}_{\mathrm{s}}|\mathrm{\Phi })\frac{\mathrm{d}{V}_{\mathrm{c}}}{\mathrm{d}z},\end{array}$$(10)

where R₀ = 20 Gpc⁻³ yr⁻¹ is the local merger rate at z = 0, ψ(z, Φ) is the merger rate evolution with redshift, and p_pop(m_s|Φ) represents the probability density function for the source frame masses m_s = (m_1, s, m_2, s). The term dV_c/dz is the differential of the comoving volume per unit redshift. For the source frame masses m_1, s and m_2, s, we adopted a power-law plus peak (PLP) model, consistent with the LIGO and Virgo GWTC-3 results (Abbott et al. 2023c):

$$\begin{array}{cc}& p(m_{1,\mathrm{s}}|m_{\min },m_{\max },\alpha )=(1-{\lambda }_{\mathrm{g}})P(m_{1,\mathrm{s}}|m_{\min },m_{\max },\alpha )\hfill \\ \hfill & +{\lambda }_{\mathrm{g}}G(m_{1,\mathrm{s}}|{\mu }_{\mathrm{g}},{\sigma }_{\mathrm{g}}),\mathrm{with}(0\le {\lambda }_{\mathrm{g}}\le 1)\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & p(m_{2,\mathrm{s}}|m_{\min },m_{1,\mathrm{s}},\beta )=P(m_{2,\mathrm{s}}|m_{\min },m_{1,\mathrm{s}},\beta ),\hfill \end{array}$$(12)

where the primary mass, m_1, s, follows a combination of a truncated power-law distribution and a Gaussian peak, weighted by the fraction λ_g. The secondary mass, m_2, s, was modeled conditionally on m_1, s using a power-law distribution. We used this model to describe a population of BBHs. The PLP model is motivated by the fact that there is a potential accumulation of CBC around 35 M_⊙ just before the pair-instability supernovae gap (Talbot & Thrane 2018). In both distributions, we applied a smoothing window at the low end of the mass distribution dependent on the parameter δ_m, which represents the smoothing length. This is introduced to model the effects of the stellar progenitor metallicity (Abbott et al. 2021c). The redshift evolution, ψ(z, Φ), was modeled using a rate model (Madau & Dickinson 2014; Callister et al. 2020), which is expressed as

$$\begin{array}{c}\hfill \psi (z,\mathrm{\Phi })=[1+{(1+z_{\mathrm{p}})}^{-\gamma -\kappa }]\frac{{(1+z)}^{\gamma }}{{(1+\frac{1+z}{1+z_{\mathrm{p}}})}^{\gamma +\kappa }},\end{array}$$(13)

where γ governs growth at low redshift, z_p marks the transition redshift, and κ determines the decline at high redshift. The hyperparameters’ values assumed for the population are reported in Table 1. According to the rate model, the total number of BBH coalescences in one year in all of the universe is about 10⁶, which agrees with existing literature (Borhanian & Sathyaprakash 2024; Ng et al. 2021).

3.1. Resolved gravitational wave sources

To simulate resolved GW sources, we generated a population of BBHs following the previously described population model. We used a simulative approach similar to Fishbach et al. (2020) and made simplified assumptions to calculate the BBHs S/N and the uncertainties for the measure of the binary parameters. For each binary, we computed the detector chirp mass (M_chirp), the symmetric mass ratio (η), and the luminosity distance (d_L). The optimal S/N (ρ_true) of the binary was calculated with the approximation used in Fishbach et al. (2020) and Mastrogiovanni et al. (2021):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{true}}={\rho }_0\theta {\left(\frac{M_{\mathrm{chirp}}}{M_{\mathrm{chirp},0}}\right)}^{5/6}\frac{d_{\mathrm{L},0}}{d_{\mathrm{L}}},\end{array}$$(14)

where ρ₀ = 12 is the reference S/N for a binary system with a chirp mass of M_chirp, 0 = 26 M_⊙ and a reference luminosity distance, d_L, 0, that is calibrated on the expected sensitivities for the O5 observing run reported in Abbott et al. (2020). For the five-detector network configuration we used d_L, 0 = 1.9 Gpc, while for the four-detector network we used d_L, 0 = 1.7 Gpc. The projection factor, θ, accounts for the binary’s orientation relative to the detector network; we defined it as the quadratic sum of the combined antenna response functions F₊ and F_× (Yunes & Siemens 2013) for both the five- and four-detector networks.

To include the effects of noise across multiple detectors, we simulated an observed S/N, ρ_det, by accounting for noise realizations in each detector. Specifically, we sampled ρ_det using a non-central χ² distribution with a non-centrality parameter, ρ_true², and 2 × N_detectors = 10 degrees of freedom. In this simulation, we applied a detection threshold of ρ_det > 12.

With the previous choices, we find that the resolvable sources correspond to about 0.05% for the five-detector configuration and 0.02% for the four-detector configuration of the overall BBH population. For this simulated scenario, removing the resolvable individual sources would decrease the SGWB by ∼6% for the five-detector network and ∼5% for the four-detector network. As a consequence, the likelihood in Eq. (1) is not separable, and that is the motivation for which our combined SGWB and CBCs scenario splits two years of GW data into two independent datasets for the resolved sources and the SGWB sources, as accounting for the resolved signals in the SGWB would require the modification of the likelihood (Zhong et al. 2025).

The resolved CBCs set contains about 1200 and 590 events detected for the five- and four-detector networks, respectively, while the set combined SGWB and CBCs contains about 600 and 295 solved events for the five- and four-detector networks, respectively. Once detections are identified, posterior samples are generated using a likelihood model described as follows:

$$\begin{array}{cc}& \mathcal{L}({\rho }_{\det },log(M_{\mathrm{chirp},\det }),{\theta }_{\det },{\eta }_{\det }|\rho ,log(M_{\mathrm{chirp}}),\theta ,\eta )\hfill \\ \hfill & =\mathcal{L}({\rho }_{\det }|\rho )\mathcal{L}(log(M_{\mathrm{chirp},\det })|log(M_{\mathrm{chirp}}))\mathcal{L}({\theta }_{\det }|\theta )\mathcal{L}({\eta }_{\det }|\eta ),\hfill \end{array}$$(15)

where ℒ(ρ_det|ρ) as we said follows a non-central χ² distribution. For M_chirp, η, and θ, we used the following normal distribution for the likelihoods:

$$\begin{array}{cc}& \mathcal{L}(log(M_{\mathrm{chirp},\det })|log(M_{\mathrm{chirp}}))=\mathcal{N}(log(M_{\mathrm{chirp}}),0.08\frac{8}{{\rho }_{\det }}),\hfill \end{array}$$(16)

$$\begin{array}{cc}& \mathcal{L}({\theta }_{\det }|\theta )=\mathcal{N}(\theta ,0.21\frac{8}{{\rho }_{\det }}),\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mathcal{L}({\eta }_{\det }|\eta )=\mathcal{N}(\eta ,0.08\frac{8}{{\rho }_{\det }}).\hfill \end{array}$$(18)

After obtaining posterior samples distributed according to this likelihood model with uniform priors, we converted them into detector frame masses and luminosity distance using the S/N conversion Eq. (14). Since there is an implied prior in this transformation, it is removed by calculating the Jacobian of the transformation from (ρ_det, M_chirp, det, η_det, θ_det) to the final variables required for our analysis (m_1, d, m_2, d, d_L).

An example of the posterior distributions for key parameters, including detector frame masses and luminosity distance, is shown in Fig. 1. With this type of likelihood, we can see that it reproduces correlated mass estimates, as expected for real GW events.

Fig. 1. Corner plot showing the posterior distributions of the key parameters estimated from the simulation (blue). The red lines represent the true values of the parameters injected into the simulation. The diagonal panels display the marginalized 1D posterior distributions for each parameter, while the off-diagonal panels show the 2D joint posterior distributions.

3.2. Stochastic gravitational-wave background

We simulated cross-correlation measurements of the corresponding SGWB, assuming T_obs = 1 yr of integration with LVK with O5 sensitivity and a coherence time for the stochastic search of 4 seconds. This corresponds to frequency point estimates of the cross-correlation statistic evaluated every 0.25 Hz. We note that this choice differs from the value of 32 seconds that is typically used for current SGWB searches. We decided to reduce the coherence time in order to speed-up our hierarchical likelihood evaluation by a factor of 8. We also note that the choice of the length of the fast Fourier transform is not expected to strongly modify our results, as the sensitivity of detection of the SGWB is almost invariant for the choice of the coherence time (Thrane & Romano 2013) and mostly dependent on the observing time. As a detector baseline we use a five-detector LIGO (Hanford, Livingston, and India), Virgo and KAGRA networks at design sensitivities¹ (Acernese et al. 2014; LIGO Scientific Collaboration 2015; Aso et al. 2013; Unnikrishnan 2013) and as a four-detector network the same but without LIGO India.

In Fig. 2, each data point is drawn from a Gaussian distribution centered on the true value of Ω_GW(f), as calculated using Eq. (4), with a standard deviation, σ(f), given by Eq. (8). The uncertainty, σ(f), represents the sensitivity to an SGWB in every frequency bin. As argued before, σ(f) is directly related to the upper limits of the energy density of the SGWB, as fluctuation beyond this value is rare. This property allows us to detect the SGWB: if σ(f) is much larger than Ω_GW(f), the variations in Ω_GW(f) become effectively indistinguishable within the noise, making detection unlikely. Conversely, when σ(f) is sufficiently small, deviations in Ω_GW(f) can be resolved, enabling us to identify the SGWB signal. Figure 2 also shows why low H₀ values are easier to exclude with the SGWB. As we can observe, when we remove the 1/H₀² dependence; namely, when we normalize the critical density of the Universe, ρ_c, the cross-correlation coefficients do not depend on H₀ (as also indicated in Eq. 5), while the SGWB does. As a result, the S/N for the SGWB increases as 1/H₀.

Fig. 2. Simulated SGWB cross-correlation spectrum, C(f) (points) with corresponding uncertainties, σ(f) (gray lines). The red lines are theoretical curves of ΩGW(f) for different h02, and the blue curves display the 1σ and 2σ sensitivity of the LVK at the design sensitivity (1 year of observation in O5). The scatter reflects realistic fluctuations centered on the true ΩGW(f).

4. Results with simulated data

Using the simulated data described in Sect. 3 and the adopted prior ranges (see Table 1), we performed three separate analyses. The first considers only resolved CBC events (resolved CBCs), using two years of data. The second is a combined analysis (combined SGWB and CBCs), where we used one year of resolved BBH detections and one year of SGWB measurements. Finally, the third analysis focuses solely on one year of the SGWB to assess its sensitivity and, if possible, directly observe the impact of improvements due to unresolved sources.

Figure 3 displays the corner plot for the cosmological parameters H₀ and Ω_m with the redshift rate parameters. For a five-detector network, from the analysis using only resolved CBCs, we find $H_0=59.7_{-13.4}^{+17.1}$ km s⁻¹ Mpc⁻¹ at 68% credible interval (CI). Combining SGWB information, we find $H_0=66.2_{-12.1}^{+19.3}$ km s⁻¹ Mpc⁻¹ at 68% CI, indicating no significant improvement in precision. This result indicates that most of the precision on H₀ is determined by resolved spectral sirens. However, we note that the inclusion of the SGWB significantly helps in excluding the region of the parameter space that covers low values of H₀ and Ω_m. As also indicated by Cousins et al. (2025), this is a consequence of the fact that the S/N for the SGWB detection is higher for low values of H₀ and Ω_m. Instead, for a four-detector network, we find that the population posterior is entirely dominated by resolved sources and the inclusion of the SGWB is not able to improve the precision in any particular region of the H₀, Ω_m plane. We assumed that this result is also valid for the rest of the population parameters. As such, for the remainder of the paper, we focus on discussing possible constraints that can be achieved by the five-detector network as we are interested in discussing the impact of the SGWB.

Fig. 3. Posterior distributions of the population parameters governing the cosmic expansion and the redshift evolution of the CBC merger rate. Results are shown for the resolved CBCs dataset (2 years of observations) using four detectors (orange) and five detectors (purple), and for the combined SGWB and CBCs dataset (1 year of individually resolved sources and 1 year of stochastic background) using four detectors (green) and five detectors (blue). Dashed red lines on the diagonal panels indicate injected (true) values.

We observe that the datasets for individual sources and SGWB measure the rate parameter, R₀, with similar precision. This is a consequence of the fact that the individual sources are mostly detected at lower redshifts and are the ones driving the inference on these population parameters. On γ we achieve an additional 10% precision with the combined analysis. We also see that the inclusion of the SGWB slightly introduces some preferences for the posteriors on z_p and k. These two parameters govern the peak of the BBH merger rate and its redshift evolution for z > z_p. The inclusion of the stochastic background makes possible the measurement $z_{\mathrm{p}}=2.8_{-0.5}^{+0.6}$ and helps to exclude low values of k; namely, we excluded models in which the rate increases beyond z_p.

Figure 4 focuses on the mass model parameters and correlations with H₀. The combined analysis does not lead to any improvement in these parameters, whereas the resolved-only analysis performs better.

Fig. 4. Posterior distributions of the population parameters governing the cosmic expansion and the CBC mass distribution. Results are shown for the resolved CBCs dataset (2 years of observations) using four detectors (orange) and five detectors (purple), and for the combined SGWB and CBCs dataset (1 year of individually resolved sources and 1 year of stochastic background) using four detectors (green) and five detectors (blue). Dashed red lines on the diagonal panels indicate injected (true) values.

To assess the performance of the inference, we produced posterior predictive checks (PPCs), which are defined as

$$\begin{array}{cc}& p(\theta |x)={\displaystyle \int p}(\mathrm{\Phi }|x)p(\theta |\mathrm{\Phi })\mathrm{d}\mathrm{\Phi },\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\mathrm{\Omega }}_{\mathrm{GW}}(f|x)={\displaystyle \int p}(\mathrm{\Phi }|x){\mathrm{\Omega }}_{\mathrm{GW}}(f;\mathrm{\Phi })\mathrm{d}\mathrm{\Phi },\hfill \end{array}$$(20)

where p(Φ|x) represents the posteriors on the population parameters and {x} the GW data that can either include or not the SGWB. By analyzing the PPCs for all three analyses in Fig. 5, we gain insight into the reconstructed redshift distribution (top panel), mass distribution (middle panel), and SGWB energy density (bottom panel), Ω_GW(f). In Fig. 5 we report the PPCs for the analyses that we have performed. When considering resolved sources alone, the reconstructed redshift distribution is precise at low redshifts, while at higher redshifts z ≥ z_p, the PPC is highly prior-dominated as we cannot determine the values of k and z_p. When using the SGWB-only data, the reconstruction of the rate is strongly prior-dominated as there is too much degeneracy among the rate parameters R₀, γ, and k, although there appears to be some information on the values of z_p. Including individual sources and SGWB together provides a reconstruction of the rate informative up to redshift z = 2 − 3 consistent with what is observed in Callister et al. (2020). As we can see, by combining the resolved sources and the SGWB, the reconstruction of the BBH merger rate is less prior-dominated.

Fig. 5. PPCs for resolved CBC analysis (left), SGWB-only analysis (center), and the joint analysis of the two sources (right). The dashed and solid blue curves show the median and 90% credible bounds, and the red curve is the fiducial value injected. Top panel: Rate density R(z) of BBH mergers over the local merger rate, R0, as a function of redshift. Middle panel: Primary mass distribution as a function of the primary mass. Bottom panel: Posterior on ΩGWh02 as a function of the frequency. In all panels, the black line represents the 1σ curve.

For the mass spectrum, we find that most of the information is collected from individual sources. In fact, for the SGWB-only analysis, the reconstruction of the mass spectrum is strongly prior dominated. The improvement in the precision of the reconstruction of SGWB energy density, when considering individual and SGWB sources, is mostly due to the improvement in the measure of γ and z_p implied by the CBC rate parameters.

5. Reanalysis of GWTC-3

Using the approach defined in Sect. 2 and the likelihood in Eq. (1), we re-performed a population and cosmological spectral siren run with data from the latest public LVK observing run. We selected 59 BBHs detected during O3 with an inverse false alarm rate higher than 1 yr and used their posteriors released with GWTC-3 (Abbott et al. 2023d). To correct selection biases for the Poisson part of the likelihood, we used a set of injections in real O3 data² released with Abbott et al. (2023c). For the stochastic data, we used the frequency point estimation of $\widehat{C}(f)$ and σ²(f) released³ with Abbott et al. (2021b). We only considered frequencies below 200 Hz to ease the computational load of the analysis. The $\widehat{C}(f)$ that we used were evaluated every 0.03 Hz, corresponding to a coherence time for the SGWB search of 32 s and they were calculated with H₀ = 67.9 km s⁻¹ Mpc⁻¹. For this analysis, we used the same BBH merger redshift model considered for the simulations. However, for the mass model we adopted a MULTI PEAK model, which can also describe a possible feature at masses of ∼10 M_⊙. The MULTI PEAK model describes the mass distributions as

$$\begin{array}{cc}& p(m_{1,\mathrm{s}}|m_{\min },m_{\max },\alpha )=(1-{\lambda }_{\mathrm{g}})P(m_{1,\mathrm{s}}|m_{\min },m_{\max },\alpha )\hfill \\ \hfill & +{\lambda }_{\mathrm{g}}{\lambda }_{\mathrm{g},\mathrm{low}}G(m_{1,\mathrm{s}}|{\mu }_{\mathrm{g},\mathrm{low}},{\sigma }_{\mathrm{g},\mathrm{low}})\hfill \\ \hfill & +{\lambda }_{\mathrm{g}}(1-{\lambda }_{\mathrm{g},\mathrm{low}})G(m_{1,\mathrm{s}}|{\mu }_{\mathrm{g},\mathrm{high}},{\sigma }_{\mathrm{g},\mathrm{high}}),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill & p(m_{2,\mathrm{s}}|m_{\min },m_{1,\mathrm{s}},\beta )=P(m_{2,\mathrm{s}}|m_{\min },m_{1,\mathrm{s}},\beta ).\hfill \end{array}$$(22)

For the reanalysis of O3, we also made the implicit assumption that the likelihood in Eq. (1) is separable. For the O3 runs, this is a reasonable assumption as the removal of the individual sources of the population does not significantly modify the SGWB (Abbott et al. 2021b). The priors used for the analysis are reported in Table 1.

We performed two analyses, one considering only individual sources and the other also considering the SGWB upper limits. We find that the constraints on the population parameters, including also the cosmological parameters, are entirely given by the individual GW sources. In Fig. 6 we display the marginalized joint posterior distributions for these two analyses on the Hubble constant, Ω_m and the parameter γ for the rate evolution. We find that the inclusion of the SGWB weakly helps in excluding low values of H₀ and high values of γ from the 2σ CI areas. The posterior on the remaining of the population parameters is unchanged and we do not display it. In Fig. 7 we show the PPCs for the BBHs rate, mass distribution and SGWB reconstructed from the two analyses. We notice that the mass and rate reconstructions are equivalent in the case that the SGWB is included or not. We can also see that the reconstructed SGWB in both analyses is well below the 1σ upper limits from the stochastic searches. This means that the SGWB prediction is entirely given by the individual sources of CBCs and no additional information can be included from the SGWB likelihood. This is consistent with what was found in Callister et al. (2020).

Fig. 6. Marginalized joint posterior distribution on the cosmological expansion parameters, and the rate parameter, γ, for an analysis including only 59 BBHs detected in O3 (purple line) and considering also the SGWB (blue line). The contour marks the 1σ and 2σ CIs.

Fig. 7. PPCs for the resolved CBC analysis (left) and the joint analysis of the two sources (right). The dashed black and solid blue curves show the median and 90% credible bounds, and the red curve is the fiducial value injected. Top panel: Rate density, R(z), of BBH mergers over the local merger rate, R0, as a function of redshift. Middle panel: Primary mass distribution as a function of the primary mass. Bottom panel: Posterior on ΩGWh02 as a function of the frequency. In all the panels, the black line represents the 1σ curve.

In parallel to this work, Cousins et al. (2025) present a spectral siren analysis using GWTC-3, including the SGWB. Although the results that we obtain are in excellent agreement, namely the inclusion of the SGWB does not modify significantly the inference on population and cosmological parameters, there are a few technical differences in our approach. From a data selection point of view, Cousins et al. (2025) used 42 BBH events with S/N > 11 observed in the first three LIGO and Virgo observations, while we used 59 BBHs detected during O3 with inverse false alarm rate higher than 1 per year. Another difference is that Cousins et al. (2025) used a PLP model to describe the masses, while we used a multi-peak model. From a data analysis perspective, Cousins et al. (2025) used as a prior on the population and cosmological parameters the posterior in Abbott et al. (2023b) from the spectral siren analysis. This was motivated by the separability of the likelihood in Eq. (1) that we also assume in our framework. Then the SGWB information is included as a likelihood term fitting posterior samples from Abbott et al. (2021b) on Ω_α = 2/3, which is a quantity describing a SGWB as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{GW}}(f)={\mathrm{\Omega }}_{\alpha =2/3}{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{2/3}\end{array}$$(23)

and corresponds to the SGWB expected by a population of BBHs emitting GWs during their inspiral phase at the 0 post-Newtonian order. Instead, we calculated Ω_GW by re-weighting a GW energy spectrum that was calculated considering corrections for the GW emission up to the 3.5 post-Newtonian order (Ajith et al. 2008) and also for the merger and ring-down phases. However, we notice that our two approaches reach equivalent results because the SGWB is dominated in the 20 − 200 Hz region by the inspiral of BBH merger that can be safely estimated with the 0 post-Newtonian order emission; therefore, in that region, the SGWB is well approximated by the prescription in Eq. (23).

6. Conclusions

In this work we have discussed how it is possible to include the SGWB in spectral siren analyses for GW cosmology and enhance our understanding of astrophysical populations and cosmological parameters. Using simulated data at the design sensitivity of the LVK network’s O5 observing run, we applied a hierarchical Bayesian framework to infer a set of 14 hyperparameters, Φ, constraining key cosmological parameters such as H₀ and Ω_m, along with parameters governing the CBC mass distribution and the redshift evolution of the merger rate.

Our analysis that uses both four- and five-detector networks shows that combining SGWB data into the inference framework does not improve the precision on the marginal posteriors for cosmological parameters. In particular, we see that the Hubble constant estimates derived solely from resolved events yield $H_0=59.7_{-13.4}^{+17.1}$ km s⁻¹ Mpc⁻¹ at 68% CI, while combining the SGWB information yields $H_0=66.2_{-12.1}^{+19.3}$ km s⁻¹ Mpc⁻¹ at 68% CI, indicating no significant improvement in precision. However, we find with the five-detector network that the inclusion of the SGWB helps exclude the lowest values of H₀ and Ω_m. Beyond cosmology, on an astrophysical population level, the integration of SGWB data slightly improves our understanding of the redshift evolution of CBC rates and helps with the measurements of γ and z_p. The SGWB’s sensitivity to unresolved sources at higher redshifts allows for tighter constraints of the merger rate distribution at high redshifts.

When analyzing real O3 data, we find that the inclusion of the SGWB does not significantly improve any of the population or cosmological parameters. This result is consistent with the findings of Callister et al. (2020) and Abbott et al. (2021c) for population studies and Cousins et al. (2025) for the cosmological expansion parameters, and it is a consequence of the fact that the SGWB implied by individual sources is well below the sensitivity estimates for current stochastic searches.

Of course, in this exploratory study, we made a few simplistic assumptions that can be revised in follow-up studies. First, to preserve the separability of the likelihood in Eq. (1), we made the suboptimal choice of splitting the data into disjoint sets used for either the SGWB or individual source searches. This can be improved by developing a likelihood that takes additional correlations with the detected sources into account. In this sense, we believe that our improvement in the determination of the cosmological parameters is “conservative” as we are not using the full information present in the dataset. We expect that, with a joint likelihood for resolved sources and the SGWB, the precision on the cosmological parameters could be improved by including the stochastic. We also did not consider the possible contribution for the SGWB introduced by binary neutron stars and neutron star black hole binaries, which can be modeled with multi-population models.

Acknowledgments

We are grateful to M. Mancarella, A. Renzini and A. Toubiana for discussions during the development of this work. We are grateful to Cousins et al. (2025) and collaborators for the constructive feedback received during the internal LVK review of this work. This work received support from the French government under the France 2030 investment plan, as part of the Excellence Initiative of Aix Marseille University – amidex (AMX-19-IET-008 – IPhU). SM is supported by ERC Starting Grant No. 101163912–GravitySirens. This research made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center, a service of the LIGO Scientific Collaboration, the KAGRA Collaboration and the Virgo Collaboration. This material is based upon work supported by NSF’s LIGO Laboratory which is a majorfacility fully funded by the National Science Foundation. The authors are grateful for computational resources provided by the LIGO Laboratory (LHO) and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

References

Appendix A: Computation of hierarchical and stochastic likelihoods

We used the Python package ICAROGW to analyze our data. For a detailed explanation of its functionality, refer to Mastrogiovanni et al. (2024). Here, we provide a summary of the likelihood computation. The main application of ICAROGW is to infer the population parameters Φ that describe the production rate of events in terms of the GW parameters, θ, namely $\frac{\mathrm{d}N}{\mathrm{d}t\mathrm{d}\theta }(\mathrm{\Phi })$. The hierarchical likelihood in Eq. 2, introduced in Sect. 2.1, models the probability of observing N_obs observations, each described by some parameters θ, in a dataset {x} over an observing time T_obs accounting for selection biases (Vitale et al. 2021). ICAROGW computes numerically the hierarchical likelihood in Eq. (2) as

$$\begin{array}{c}\hfill ln[\mathcal{L}(\{x\}|\mathrm{\Phi })]\approx -\frac{T_{\mathrm{obs}}}{N_{\mathrm{gen}}}{\displaystyle \sum _{j=1}^{N_{\det }}}{s}_j+{\displaystyle \sum _{i=1}^{N_{\mathrm{obs}}}}ln\left[\frac{T_{\mathrm{obs}}}{N_{s,i}}{\displaystyle \sum _{j=1}^{N_{s,i}}}{w}_{i,j}\right],\end{array}$$(A.1)

where w_i, j and s_j are weights: w_i, j has a dimension equal to the number of CBC mergers happening per unit of time and weights the sample, while s_j weights the injections and it is defined with the dimension of a CBC merger rate per detector time. The first term derives from Monte Carlo integration to estimate

$$\begin{array}{c}\hfill N_{\exp }(\mathrm{\Phi })=T_{\mathrm{obs}}{\displaystyle \int p_{\det }}(\theta ,\mathrm{\Phi })\frac{\mathrm{d}N}{\mathrm{d}t\mathrm{d}\theta }\mathrm{d}\theta ,\end{array}$$(A.2)

where p_det(θ, Φ) is a detection probability or selection bias. Since an analytical form of p_det(θ, Φ) is typically unavailable, selection biases are evaluated using Monte Carlo simulations of injected and detected events (Abbott et al. 2021c). ICAROGW takes in input a set of N_det detected injections out of N_gen total injections generated from a prior π_inj(θ) to calculate the selection bias.

For each population model, ICAROGW calculates two numerical stability estimators. The first is the effective number of posterior samples per event:

$$\begin{array}{c}\hfill N_{\mathrm{eff},i}=\frac{{\left({\displaystyle {\sum }_{j=1}^{N_{s,i}}}{w}_{i,j}\right)}^2}{{\displaystyle {\sum }_{j=1}^{N_{s,j}}}{w}_{i,j}^2},\end{array}$$(A.3)

which ensures sufficient samples contribute to the integral, with a typical threshold of at least 20 effective samples per event (Callister et al. 2022; Talbot & Golomb 2023). The second is the effective number of injections:

$$\begin{array}{c}\hfill N_{\mathrm{eff},\mathrm{inj}}=\frac{{\left[{\displaystyle {\sum }_{j=1}^{N_{\det }}}{s}_j\right]}^2}{[{\displaystyle {\sum }_{j=1}^{N_{\det }}}{s}_j^2-N_{\mathrm{gen}}^{-1}{\left(\left[{\displaystyle {\sum }_{j=1}^{N_{\det }}}{s}_j\right]\right)}^2]},\end{array}$$(A.4)

where numerical stability typically requires N_eff, inj > 4N_obs (Callister et al. 2022; Talbot & Golomb 2023).

The likelihood computation involves technical flags for configuration: nparallel = 2048, which is the number of posterior samples that will be used per event to compute w_i, j, neffPE = −1, which indicates that there is no threshold for the effective number of posterior samples per event in Eq. A.3 and neffINJ=None, which is the effective number of injections required by the hierarchical likelihood. If neffINJ=None, a default threshold of NeffINJ = 4 × N_obs will be used.

We developed and implemented in ICAROGW a new function to calculate the SGWB likelihood, as described in Eq. 9, optimizing the computation of Ω_GW. The expected Ω_GW(f) calculation, conditioned on a set of hyperparameters, Φ, follows the form outlined in Eq. 4. This value depends on the average energy, $\langle \frac{\mathrm{d}E}{\mathrm{d}{f}_s}\rangle$, radiated by individual binaries, which is expressed as

$$\begin{array}{c}\hfill \langle \frac{\mathrm{d}E}{\mathrm{d}{f}_s}\rangle ={\displaystyle \int \mathrm{d}}\mathrm{\Phi }p(\mathrm{\Phi })\frac{\mathrm{d}{E}_s}{\mathrm{d}{f}_s}(\mathrm{\Phi }),\end{array}$$(A.5)

where Φ denotes the intrinsic properties of a given binary (masses, redshift, etc.), p(Φ) is the distribution of these parameters across the CBC population, and this quantity is evaluated at the source-frame frequency f_s = f(1 + z) in Eq. 4. The GW energy spectrum is calculated with functionalities from PYGWB Renzini et al. (2023) and Turbang et al. (2024) that calculated the GW energy spectrum up to the 3.5 post-Newtonian order for the GW emission and the merger and ring-down phase for non-precessing binaries.

Here, the integration spans the component masses m₁ and m₂. Extending this to include integration over redshift, the computation of the expected Ω_GW(f) becomes a 3D integral. A practical alternative is to use Monte Carlo integration following the methodology used in Turbang et al. (2024), which estimates the integral by averaging over numerous samples drawn from the population.

In the implementation of this approach, we started by generating N random samples from uniform distributions for the parameters z, m₁, and m₂. Let the i-th sampled set of parameters be denoted as {z_i, m_1, i, m_2, i}. For each of these samples, we calculated the energy spectrum, denoted as ${\left(\frac{\mathrm{d}E}{\mathrm{d}{f}_s}\right)}_i$. However, to evaluate Ω_GW(f) for the population described by the hyperparameters Φ, re-weighting of these samples is necessary. The computation is expressed schematically as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{GW}}(f)\propto \frac{1}{\mathrm{N}}{\displaystyle \sum _{i=1}^N}{w}_i{\left(\frac{\mathrm{d}E}{\mathrm{d}{f}_s}\right)}_i,\end{array}$$(A.6)

where w_i are weights defined as

$$\begin{array}{c}\hfill w_i\propto \frac{\mathcal{R}(z_i)(1+z_i)H{(z_i)}^{-1}p(m_{1i})p(m_{2i})}{p_{\mathrm{draw}}(z_i)p_{\mathrm{draw}}(m_{1i})p_{\mathrm{draw}}(m_{2i})}.\end{array}$$(A.7)

Here, p_draw refers to the uniform distributions from which the initial samples were drawn, and the numerators represent the target distributions associated with the population characterized by Φ. This weighting allows us to compute Ω_GW(f) while accounting for the desired population properties.

Appendix B: Measurement of the expansion rate, H(z)

To better understand the impact of combining resolved and unresolved sources on the reconstruction of the cosmic expansion history. We observe a modest improvement in the measurement of the expansion rate, H(z), when including the SGWB in the five-detector analysis, as shown on the left in Fig. B.1. This improvement is driven by the additional high-redshift information provided by unresolved sources, which is not accessible through individually resolved events alone. In particular, the SGWB contribution enables us to better constrain the high-redshift behavior of the expansion by disfavoring low values of both Ω_m and H₀, leading to a slightly tighter reconstruction of the expansion history. In contrast, the four-detector analysis, shown on the right, reveals no appreciable difference between the resolved-only and combineddatasets.

Fig. B.1. Posterior predictive distributions for the resolved CBCs dataset (purple) and the combined SGWB and CBCs dataset (blue) for results based on five (right) and four detectors (left). The solid black lines indicate the median predictions, with dashed lines showing the 90% CIs. The red curve denotes the injected (fiducial) model.

All Tables

All Figures

	Fig. 1. Corner plot showing the posterior distributions of the key parameters estimated from the simulation (blue). The red lines represent the true values of the parameters injected into the simulation. The diagonal panels display the marginalized 1D posterior distributions for each parameter, while the off-diagonal panels show the 2D joint posterior distributions.
In the text

	Fig. 2. Simulated SGWB cross-correlation spectrum, C(f) (points) with corresponding uncertainties, σ(f) (gray lines). The red lines are theoretical curves of ΩGW(f) for different h02, and the blue curves display the 1σ and 2σ sensitivity of the LVK at the design sensitivity (1 year of observation in O5). The scatter reflects realistic fluctuations centered on the true ΩGW(f).
In the text

Fig. 3. Posterior distributions of the population parameters governing the cosmic expansion and the redshift evolution of the CBC merger rate. Results are shown for the resolved CBCs dataset (2 years of observations) using four detectors (orange) and five detectors (purple), and for the combined SGWB and CBCs dataset (1 year of individually resolved sources and 1 year of stochastic background) using four detectors (green) and five detectors (blue). Dashed red lines on the diagonal panels indicate injected (true) values.

In the text

Fig. 4. Posterior distributions of the population parameters governing the cosmic expansion and the CBC mass distribution. Results are shown for the resolved CBCs dataset (2 years of observations) using four detectors (orange) and five detectors (purple), and for the combined SGWB and CBCs dataset (1 year of individually resolved sources and 1 year of stochastic background) using four detectors (green) and five detectors (blue). Dashed red lines on the diagonal panels indicate injected (true) values.

In the text

Fig. 5. PPCs for resolved CBC analysis (left), SGWB-only analysis (center), and the joint analysis of the two sources (right). The dashed and solid blue curves show the median and 90% credible bounds, and the red curve is the fiducial value injected. Top panel: Rate density R(z) of BBH mergers over the local merger rate, R0, as a function of redshift. Middle panel: Primary mass distribution as a function of the primary mass. Bottom panel: Posterior on ΩGWh02 as a function of the frequency. In all panels, the black line represents the 1σ curve.

In the text

	Fig. 6. Marginalized joint posterior distribution on the cosmological expansion parameters, and the rate parameter, γ, for an analysis including only 59 BBHs detected in O3 (purple line) and considering also the SGWB (blue line). The contour marks the 1σ and 2σ CIs.
In the text

Fig. 7. PPCs for the resolved CBC analysis (left) and the joint analysis of the two sources (right). The dashed black and solid blue curves show the median and 90% credible bounds, and the red curve is the fiducial value injected. Top panel: Rate density, R(z), of BBH mergers over the local merger rate, R0, as a function of redshift. Middle panel: Primary mass distribution as a function of the primary mass. Bottom panel: Posterior on ΩGWh02 as a function of the frequency. In all the panels, the black line represents the 1σ curve.

In the text

	Fig. B.1. Posterior predictive distributions for the resolved CBCs dataset (purple) and the combined SGWB and CBCs dataset (blue) for results based on five (right) and four detectors (left). The solid black lines indicate the median predictions, with dashed lines showing the 90% CIs. The red curve denotes the injected (fiducial) model.
In the text