The Metallicity Dependence and Evolutionary Times of Merging Binary Black Holes: Combined Constraints from Individual Gravitational-wave Detections and the Stochastic Background

The advent of gravitational-wave astronomy is now allowing for the study of compact binary merger demographics throughout the Universe. This information can be leveraged as tools for understanding massive stars, their environments, and their evolution. One active question is the nature of compact binary formation: the environmental and chemical conditions required for black hole birth and the time delays experienced by binaries before they merge. Gravitational-wave events detected today, however, primarily occur at low or moderate redshifts due to current interferometer sensitivity, therefore limiting our ability to probe the high-redshift behavior of these quantities. In this work, we circumvent this limitation by using an additional source of information: observational limits on the gravitational-wave background from unresolved binaries in the distant Universe. Using current gravitational-wave data from the first three observing runs of LIGO–Virgo–KAGRA, we combine catalogs of directly detected binaries and limits on the stochastic background to constrain the time-delay distribution and metallicity dependence of binary black hole evolution. Looking to the future, we also explore how these constraints will be improved at the Advanced LIGO A+ sensitivity. We conclude that, although binary black hole formation cannot be strongly constrained with today’s data, the future detection (or a nondetection) of the gravitational-wave background with Advanced LIGO A+ will carry strong implications for the evolution of binary black holes.


INTRODUCTION
Since the first detection of gravitational waves by the LIGO-Virgo collaborations in 2015 (Aso et al. 2013;Aasi et al. 2015;Acernese et al. 2015;Abbott et al. 2016), the number of detected compact binary coalescences has been steadily increasing, amounting to around 90 detections in the third observing run (Abbott et al. 2023a).These individual detections offer new insights into the demographics of compact binary mergers, including their mass and spin distributions.At the same time, growing gravitational-wave catalogs offer an opportunity to study the evolution of the binary merger rate over cosmic time.The redshift-dependent merger rate is dictated by a confluence of environmental conditions governing compact binary formation and evolutionary processes describing their subsequent evolution.It therefore offers a window through which we can learn kevin.turbang@vub.beabout the births of compact binaries and the lives of their progenitors.
Compact binaries are expected to preferentially arise from massive stars in low metallicity environments (Fryer et al. 2012;Belczynski et al. 2016;Chruslinska et al. 2018;Mapelli et al. 2019;Santoliquido et al. 2020;Liotine et al. 2023).Both the underlying metallicityspecific star formation rate and the metallicity dependence of compact object formation, however, remain poorly understood (Chruslinska & Nelemans 2019;Broekgaarden et al. 2021;Broekgaarden et al. 2022;Chruślińska 2022;Chruślińska et al. 2023;van Son et al. 2023).Therefore, gravitational-wave observations might provide one of the best observational routes to constraining or measuring the metallicity specific star formation rate and compact object formation efficiency.Furthermore, the time-delay distribution between progenitor formation and eventual binary merger is not known and, if measured, it could provide hints as to what are the dominant formation channels and environments of binary black holes.Binaries in the field, for example, exhibit longer/shorter evolutionary time delays if they undergo stable/unstable (e.g. common envelope) mass transfer, respectively (van Son et al. 2022).Binaries assembled dynamically in clusters might experience a different range of time delays altogether, with further distinctions depending on whether these dynamical systems merge inside the cluster or are jettisoned and merge outside the cluster (Samsing & D'Orazio 2018).
Current gravitational-wave detectors, however, are limited in their ability to probe high-redshift binary black hole merger events.For example, current catalogs provide meaningful constraints on the merger rate only to z ≲ 1 (Abbott et al. 2023a;Callister & Farr 2023), although the formation rate can be probed out to z ∼ 4 for population synthesis-motivated assumptions about the time-delay distribution (Fishbach & van Son 2023).Even with the future Advanced LIGO A+ sensitivity (Barsotti et al. 2018), the redshifts that will be probed by individual compact binary detections will likely remain below z ∼ 3, thus suggesting that the highredshift investigation of the binary black hole merger rate will remain challenging, at least until the era of next-generation instruments like the Einstein Telescope and Cosmic Explorer (Vitale et al. 2019;Evans et al. 2021;Borhanian & Sathyaprakash 2022).
Direct compact binary detections are not our only source of information, however.In addition to those individual binaries detected by the LIGO-Virgo-KAGRA collaborations, a gravitational-wave background from unresolvable compact binaries is expected (Romano & Cornish 2017;Christensen 2018;van Remortel et al. 2023).This gravitational-wave background can be detected in the form of excess correlated noise shared between distant gravitational-wave detectors.Even though a gravitational-wave background remains to be detected by the LIGO-Virgo-KAGRA collaborations, current upper limits on the gravitational-wave background can already allow us to constrain the merger rate at higher redshifts than individual binary black hole merger detections alone (Abbott et al. 2021a,b).
Many studies have previously used catalogs of gravitational-wave detections to study the formation history of compact binaries (Fishbach et al. 2018;Vitale et al. 2019;Fishbach & Kalogera 2021;Mukherjee & Dizgah 2022;Riley & Mandel 2023).In this paper, we will study the evolutionary history of compact binaries by combining direct detections with constraints on the gravitational-wave background.The approach of combining individual binary black hole merger events with upper limits on the gravitational-wave background was first explored in Callister et al. (2020), where a phenomenological broken power-law model was adopted for the binary black hole merger rate.In this study, how-ever, we will seek to go a step further, moving past the merger rate itself and constraining the metallicity dependence and evolutionary time delays associated with compact binary formation and evolution.
We find that current data can constrain the slope of the binary time-delay distribution, but does not contain information regarding the metallicity dependence of compact binary formation.Moreover, current limits on the stochastic gravitational-wave background are not yet strong enough to noticeably impact these results.However, we show that the sensitivity of future A+ instruments will dramatically change this picture.At A+ sensitivities, both the detection or non-detection of the gravitational-wave background will provide novel constraints on the parameters governing these distributions, offering a unique probe of the environment in which the binary formed and its evolution.
This work is structured as follows.In Sec. 2, we introduce our model for the binary black hole merger rate, starting from a metallicity-dependent star formation rate and an evolutionary time-delay distribution.In parallel, we demonstrate how differing assumptions regarding binary black hole formation impact the expected amplitude of the gravitational-wave background.The methodology for our joint individual binary black hole and stochastic background analysis method is outlined in Sec. 3. Our results using data from the first three LIGO-Virgo-KAGRA observing runs are discussed in Sec.4.1.Finally, we explore the potential of the analysis method with future constraints by considering the Advanced LIGO A+ sensitivity in Sec.4.2.

BINARY BLACK HOLE MERGER RATE
Whereas past work (Callister et al. 2020) focused on characterizing the resulting merger rate of black holes, in this paper we will seek to work further backwards, studying the underlying formation efficiency and timedelay distribution of binary black holes.Standard models for the merger history of compact binaries rely on three ingredients: a metallicity specific star formation rate, an efficiency with which stellar progenitors yield merging compact objects, and a distribution of evolutionary time delays.
Throughout this work, we will consider the star formation rate as given by Madau & Dickinson (2014).They obtain a fit to the global star formation rate per unit volume of the form: 1 + ((1 + z)/2.9)5.6 . (1) In App.A, we provide results for an alternative star formation rate given by Vangioni et al. (2015) to illus-trate that the main conclusions of this work can still be drawn if one assumes a different star formation rate.
It is expected that merging binary black holes preferentially form in low metallicity environments, with efficiencies falling steeply at metallicities Z ≳ 0.1 Z ⊙ (Belczynski et al. 2016;Chruslinska et al. 2018;Klencki et al. 2018;Mapelli et al. 2019;Neijssel et al. 2019;Santoliquido et al. 2020;Santoliquido et al. 2021;Broekgaarden et al. 2022).Motivated by this, we will assume that the birth rate of binary black holes is proportional to the product R * (z)F (Z max , z), where F (Z max , z) is the fraction of star formation occurring at redshift z below metallicity Z max .An analytic prescription for the fraction of star formation is given by Langer & Norman (2006) as where Γ is the incomplete gamma function and Z ⊙ is the solar metallicity.The merger rate as a function of redshift is then obtained by performing a convolution of the metallicity-dependent star formation rate and a distribution of time delays t d between binary formation and merger.We model this time-delay distribution as a power law, where t max d is fixed to the age of the Universe, i.e., 13.5 Gyrs, and t min d is a parameter that can be inferred from the data.The merger rate as a function of redshift is now given by the convolution of the time-delay distribution with the metallicity-dependent star formation rate: where the binary formation redshifts z f = z f (z, t d ) are regarded as a function of merger redshift and time delay.
Our prescription for the binary black hole merger rate left three parameters undefined -the slope κ of the timedelay distribution, the minimum time delay t min d , and the maximum metallicity Z max amenable to compact binary formation.Our goal in this work will be to infer these parameters from gravitational-wave observations.The measurement of these parameters could shed light on the different binary formation channels, as each of these predicts different values of the spectral index κ.For example, the value of κ is expected to be -1 in the classical field formation scenario (Dominik et al. 2012;Dominik et al. 2013;Fishbach & van Son 2023).
An example of the merger rate R(z) is given in the left-hand side of Figure 1.We vary the slope of the timedelay distribution κ, as well as the maximum metallicity below which black hole formation occurs Z max , but fix the value of the minimum time-delay parameter t min d to an arbitrary value of 0.05 Gyrs.Each curve is arbitrarily normalized to R(z = 0.2) = 1.A more negative slope κ corresponds to more events that merged at early evolutionary times.Given a fixed observed merger rate in the local Universe, more negative κ correspondingly boosts the merger rate at large redshifts, as illustrated in the left-hand side of Figure 1.Similarly, lower metallicity thresholds allow for a larger merger rate at earlier times, i.e., at large redshifts, when the metallicity of the environment in which the binaries formed was still low.Note that, although the effect of varying the minimum time-delay parameter is not displayed, decreasing its value would result in an increased merger rate at high redshifts, as mergers now occur at earlier times, and therefore, larger redshifts.
As we vary parameters, the behavior of the black hole merger history in the local Universe, where binary black hole mergers are currently observed, remains practically unaltered (Abbott et al. 2023a).In contrast, the most dramatic differences occur at large redshifts.Although these redshifts are too distant to directly observe, these different models yield very different predictions regarding the stochastic gravitational-wave background.
The gravitational-wave background is expressed in terms of the dimensionless energy density fraction (Romano & Cornish 2017): where ρ c = 3H 2 0 c 2 /(8πG) is the critical energy density of the Universe, G is Newton's gravitational constant, c is the speed of light, and the Hubble constant is given by H 0 = 67.9km s −1 Mpc −1 (Planck Collaboration et al. 2016).More specifically, for the gravitational-wave background coming from unresolved compact binaries throughout the Universe, the dimensionless energy density takes the form (Phinney 2001;Callister et al. 2020): where the Hubble rate at redshift z is given by Redshift z The value of all other parameters that enter in the mass distribution as given by Eqs. ( 14) -( 16) is fixed to arbitrary values: (3) and ( 12), respectively.The behavior of the Ω(f ) spectrum as the time-delay parameters are varied can directly be related to the behavior of the merger rate on the left-hand side, as a larger integrated merger rate results in a larger net Ω(f ) spectrum.
holes are expected beyond this redshift, provided the binary black holes are stellar progenitors.The quantity ⟨dE s /df s ⟩ is the gravitational-wave energy spectrum of a single binary, averaged across the binary black hole population (Ajith et al. 2008).If λ denotes the intrinsic properties of a given binary (masses, spins, etc.) and p(λ) is the distribution of these parameters across the binary black hole population, then where this quantity is evaluated at the source-frame frequency f s = f (1 + z) in Eq. ( 6).More details about the assumed mass and spin distributions will be given below.
In the right-hand side of Figure 1, we show an example of the computed Ω(f ) spectrum for different values of the time-delay distribution parameters, as introduced in Eqs. ( 3) and ( 4).The variation in the Ω(f ) spectrum when varying the time-delay distribution and metallicity parameters can be directly linked to the behavior of the merger rate R(z) in the left-hand side of Figure 1.Indeed, increasing the integrated merger rate correspondingly increases the Ω(f ).Note that the largest background comes from a steep time-delay distribution and a very low maximum metallicity, while the smallest background comes from allowing long time delays and a high maximum metallicity.For a more detailed expla-nation of the computation of Ω(f ), we refer the reader to App.D. Furthermore, we note that the example given here is largely illustrative.More detailed investigations have recently been performed by Lehoucq et al. (2023).

ANALYSIS METHOD
Our goal is to learn about the evolutionary history of binary black holes by synthesizing all available gravitational-wave information, combining both direct detections of black hole mergers with upper limits on the integrated gravitational-wave background, based on work by Callister et al. (2020).The inputs to our analysis are a set of N obs direct binary black hole detections, with data {d} N obs i=1 , along with cross-correlation measurements Ĉ(f ) of the gravitational-wave background (to be explained further below).Let Λ signify the set of hyperparameters describing the binary black hole population, e.g., the slope of the time-delay distribution, κ, introduced in Eq. (3).We assume that the joint likelihood of these measurements can be factorized as where p BBH and p GWB denote the likelihoods for the individual binary black hole detections and the gravitational-wave background, respectively. 1 The individual binary black hole likelihood takes the form (Loredo 2004;Taylor & Gerosa 2018;Mandel et al. 2019): where, as in Sec. 2, λ denotes the parameters (redshift, masses, spins, etc.), and p(d i |λ) denotes the likelihood of a gravitational-wave event d i given parameters λ.In the expression above, we also note the presence of the differential mass-redshift distribution for binary black holes, which is given by: where N denotes the total number of binary black hole mergers integrated across all redshifts and masses.Assuming that the component mass distributions do not vary with redshift, 2 this allows us to decompose the mass-redshift distribution and adopt the following parameterization: where is the differential source-frame merger rate at redshift z = 0.2 and primary component mass m 1 = 20M ⊙ , and ϕ(m 1 ) and p(m 2 ) characterize the primary and secondary mass distributions, respectively.The sourceframe merger rate in Eq. ( 11) can be related to the detector-frame rate in Eq. ( 10) by where the factor (1 + z) −1 transforms from detectorframe to source-frame time, and dVc dz denotes the comoving volume per unit redshift.
1 This factorization does not strictly hold, as the same stretches of data contribute to both the direct compact binary detections and the gravitational-wave background measurements.However, the contributions of direct detections to gravitational-wave background measurements are currently negligible (Abbott et al. 2021a), making the factorization a good approximation. 2 Although intrinsic evolution of the compact binary mass distribution with redshift is a generic prediction (van Son et al. 2022;Ye & Fishbach 2024), no such evolution is detected with current data (Fishbach et al. 2021;Abbott et al. 2023a;van Son et al. 2022).
Following results from the LIGO-Virgo-KAGRA GWTC-3 catalog (Abbott et al. 2023a), we model the primary mass distribution as a mixture between a power-law and an additional Gaussian excess: where α is the spectral index of the power-law, the mean µ m and variance σ 2 m characterize the Gaussian peak, and f p denotes the relative contribution of the Gaussian peak and the power-law (Talbot & Thrane 2018;Abbott et al. 2021c).A smoothing function is applied, such that the mass distribution is suppressed for m 1 < m low and m 1 > m high : where δm low and δm high denote the scale of smoothing.The secondary mass m 2 is assumed to follow a powerlaw distribution, while ensuring m 2 < m 1 : where β q denotes the spectral index of the distribution, and p(m 2 |m 1 ) is set to zero below 2M ⊙ .Although not explicitly mentioned throughout this section, we additionally measure the distributions of component spin magnitudes and spin-orbit tilt angles; see App.B. In Eq. ( 9), the total number of gravitational-wave events, i.e., detected and undetected, is denoted by N (Λ), whereas the number of observed events is represented by N obs .The expected number of observations N exp is given by where P det (λ) is the probability of successfully detecting a binary black hole merger described by parameters λ.In practice, N exp (Λ) can be approximated through a Monte Carlo average over artificial signals injected into detector data: where the sum runs over the number of above detection threshold injections, N found , drawn from a reference distribution p inj (λ i ) (Callister et al. 2022).More details about the set of injections used in this analysis can be found in App. C. The remaining necessary element is the integration over individual binary black hole likelihoods p(d i |λ) in Eq. ( 9).This can be estimated using Monte Carlo integration over posterior samples drawn from a reference prior (Loredo 2004;Taylor & Gerosa 2018;Mandel et al. 2019): where p pe (λ i ) denotes the prior used during the parameter estimation to infer the distribution of the parameters describing the gravitational-wave event, and the average is taken over the posterior samples of each event.
Our joint analysis also incorporates the contribution from the gravitational-wave background.Gravitationalwave background searches, such as the one performed by the LIGO-Virgo-KAGRA collaborations (Abbott et al. 2021a), aim to measure the dimensionless energy density Ω(f ) introduced above.To this end, an optimal crosscorrelation estimator is defined (Allen & Romano 1999): where T is the observation time, and s1 (f ) and s2 (f ) denote the Fourier transformed data in both detectors.
The overlap reduction function γ 12 (f ) encodes the baseline geometry of the detector pair (Christensen 1992;Flanagan 1993).The normalization of this estimator is such that ⟨ Ĉ(f )⟩ = Ω(f ).An estimator of the variance is given by Ĉ under the assumption of a small signal-to-noise ratio, where P 1 (f ) and P 2 (f ) denote the one-sided power spectra in both detectors, and ∆f is the frequency resolution.This cross-correlation estimator then enters in the likelihood p GWB in Eq. ( 8), which is well-approximated by a Gaussian distribution (Abbott et al. 2021a): where the sum runs over discrete frequency bins f k and Ω(f ) is given by Eq. ( 6).For additional information, we refer the reader to Abbott et al. (2021a).
All together, we infer 16 hyperparameters: the nine parameters governing the primary and secondary mass distributions in Eqs. ( 15) and ( 16), the three spin parameters detailed in App.B, the three parameters {Z max , κ, t min d } characterizing the formation and evolution of binary black hole systems, as well as the reference merger rate amplitude R ref defined in Eq. ( 12).The priors used for the parameter estimation are summarized in Table 1 in App.B.

RESULTS
In this section, we discuss the results of our analyses.Sec.4.1 presents the results using currently available LIGO-Virgo data, while Sec.4.2 shows results that will be available in future observing runs with the improved A+ sensitivity of the Advanced LIGO detectors.

Constraints from LIGO-Virgo O1, O2, and O3 data
We now apply the methodology outlined in the previous sections to place constraints on the time-delay parameters giving rise to the binary black hole population, using data from LIGO-Virgo's first three observing runs.More details about the exact data used are given in App. C.
In Figure 2, we show the obtained posteriors on the time-delay parameters {κ, Z max , t min d } from Eqs. (3) and (4), as well as the merger rate amplitude R ref at z = 0.2 and m 1 = 20M ⊙ , as defined by Eq. ( 12).We simultaneously infer all the hyperparameters, e.g. also appearing in the mass distributions, but these are omitted from the plot for the sake of visual clarity.Since the merger rate amplitude R ref parameterizes the merger rate at z = 0.2 and m 1 = 20M ⊙ , it will be almost solely determined by LIGO-Virgo's individual binary black hole detections, which happen at low to moderate redshifts.This is indeed the case, as illustrated by the recovery of the R ref parameter in Figure 2, which is dominated by individual binary black hole events.In addition, the slope of the time-delay distribution κ is constrained to negative values.A positive value of κ would give rise to a merger rate R(z) that decreases with z at low redshifts, contradicting the currently observed merger rate dependence on the redshift.Current individual binary black hole detections together with the gravitational-wave background estimator are not able to constrain the metallicity threshold Z max , nor the minimum time-delay parameter t min d .Fishbach & Kalogera (2021) similarly constrained the conditions of compact binaries using direct detections of binary mergers from the first two LIGO observing runs (but without additional constraints from the gravitational-wave background).We do note, however, that they assume an alternate star formation rate given by Madau & Fragos (2017), and that only two parameters are varied at a time, contrarily to our case where the slope, the minimum time delay, and the metallicity threshold are varied simultaneously.Keeping these differences in mind, our results are consistent with the ones obtained by Fishbach & Kalogera (2021).They constrained the slope of the time delay distribution to negative values at the 90% confidence level, and are unable to constrain the value of the metallicity threshold for the values of the minimum time delay considered in this paper.
In principle, measurements or upper limits on the gravitational-wave background are dominated by binaries at higher redshift, and could thus offer a complementary constraint on the slope of the time-delay distribution κ, the metallicity threshold Z max , and the minimum Posterior samples on the merger rate dR dm 1 (m1 = 20M⊙) for the joint analysis using both direct detections and stochastic background constraints.For reference, the 5% and 95% credible limits are shown in full black lines for the individual binary black hole run, and the dashed black line shows the median value for the binary black hole only run.The results of the individual binary black hole run and the joint analysis using both direct detections and stochastic background constraints are comparable, illustrating that the gravitational-wave background upper limits are not informative at current detector sensitivity.Furthermore, we note that the merger rate is well-constrained at low redshifts because of the individual binary black hole mergers detected at those redshifts.Right: Posterior on Ω(f ) using the posterior samples from the joint analysis using both direct detections and stochastic background constraints on the first three LIGO-Virgo-KAGRA observing runs.The median is denoted by the dashed black line, and the 95% credible upper limit by the full black line.The 2σ O3 PI curve from Abbott et al. ( 2021a) is shown (dotted black) as an indication of the gravitational-wave background detection threshold after the first three LIGO-Virgo observing runs.All posterior samples fall below the O3 PI curve, illustrating that we do not expect the gravitational-wave background contribution to be informative at current detector sensitivity.time delay t min d .However, we find that current constraints on the gravitational-wave background are not sufficiently sensitive to appreciably change the measurements using solely individual binary black holes.To illustrate this point further, we compare the resulting posterior samples on the merger rate dR dm1 (m 1 = 20M ⊙ ) in the left-hand side of Figure 3.The individual blue traces correspond to individual posterior samples drawn from Figure 2, while solid black curves show the bounds obtained from analyzing binary black holes only.This shows that, at current sensitivity, results are entirely dominated by individual binary black hole mergers.
To illustrate why current upper limits on the gravitational-wave background are not yet informative, we show the expected Ω(f ) from the posterior samples of the joint analysis using both direct detections and stochastic background constraints in the right-hand side of Figure 3.We compare these samples with the 2σ O3 power-law integrated sensitivity curve (PI curve), which illustrates the sensitivity to an isotropic gravitationalwave background (Thrane & Romano 2013).Approximately, spectra lying above the PI curve are typically detectable with a signal-to-noise ratio ≥ 1, while spectra lying below have a signal-to-noise ratio < 1.As can be seen, none of the posterior samples for Ω(f ) are within reach of the PI curve, illustrating that we do not expect additional information from the gravitational-wave background contribution to the likelihood.
Note that our results rely on our fairly strong systematic choice of a fixed star formation rate.As a check on the corresponding systematic uncertainty inherent in our results, we repeat this analysis in App.A for an alternative star formation rate given by Vangioni et al. (2015), illustrating that our main conclusions remain unaltered when assuming a different star formation rate.The stellar metallicity distribution in Eq. ( 2) represents an additional systematic uncertainty in our analysis.Although we use a fit provided by Langer & Norman (2006) in this work, more modern fits have recently been proposed (e.g., Chruślińska et al. 2021;van Son et al. 2023).Future work will assess the impact of the choice of metallicity distribution on our results.

Constraints using Advanced LIGO A+ sensitivity
Although current gravitational-wave background limits are not yet informative, increased integration time will give increasingly sensitive gravitational-wave background measurements, even in cases where detector horizons do not appreciably change.This increased sensitivity will be increasingly informative regarding the evolutionary history of binary black holes.We therefore consider stochastic background measurements obtain-able with the Advanced LIGO detectors at their future A+ sensitivities (Barsotti et al. 2018), assessing the potential of our joint analysis method in constraining the time-delay parameters.
Following a future observing run with A+ LIGO instruments, there are two possibilities: a detection of the gravitational-wave background could be made, or the background could go undetected still.We will consider both cases, exploring the resulting constraints on the time-delay parameter space in each case.In both cases, we simulate a gravitational-wave background signal consistent with the results obtained in Sec.4.1 using data from the first three observing runs.To accomplish this, we choose two samples from the posterior distribution presented in Figure 2: one posterior sample predicting a stochastic background that should be detectable with A+ LIGO, and a second predicting a background that should remain undetected.For both cases, we consider the LIGO Hanford and LIGO Livingston baseline, and assume an observation time of one year at the LIGO Advanced A+ sensitivity (Barsotti et al. 2018).Future observing runs with more sensitive instruments will certainly also yield additional direct detections of binary black hole mergers.Our goal in this section, however, is primarily to understand the astrophysical information contained in a detection (or non-detection) of the stochastic background.We therefore do not simulate additional compact binary mergers, but use the same GWTC-3 catalog of direct detections as in Sec.4.1.
Case 1: A detection -We first consider the detectable gravitational-wave background case, for which the simulated Ω(f ) spectrum is shown in Figure 4 (right), together with the 2σ A+ PI curve assuming a oneyear observation time, for reference.As the injected gravitational-wave background lies above the PI curve, we indeed expect a detection.We infer the slope of the time-delay distribution κ, the metallicity threshold Z max , and minimum time delay t min d using both the individual binary black hole events from the GWTC-3 catalog, and the simulated detectable gravitationalwave background at A+ sensitivity.The posteriors for these parameters are shown in Figure 5, where the black dashed histograms represent the results of the analysis with current data reported in Figure 2 to facilitate comparison.We omit the R ref posterior from the plot, and note that since the merger rate amplitude R ref is currently entirely determined by individual binary black hole events, the addition of gravitational-wave background measurement at A+ sensitivity does not alter the posterior on R ref compared to the runs at current sensitivity in the previous section.Contrarily, the posteriors on the time-delay parameters are much more con- for the joint analysis using both direct detections and stochastic background constraints at Advanced LIGO A+ sensitivity, in the presence of a detectable gravitational-wave background.The dashed black lines represent the 1D histograms for the joint analysis on data from the first three observing runs from Figure 2, for reference.We note that there is more support for larger negative values of the slope of the time-delay distribution κ, and smaller values of the metallicity parameter Z max and the minimum time delay t min d are favored.In addition, we point out the complementarity of these constraints with the case of an undetectable gravitational-wave background at A+ sensitivity, as reported in Figure 7.
strained than in the runs using data from the first three observing runs of the previous section.Larger negative values are favored for the slope of the time-delay distribution κ, whereas smaller values are supported for both the metallicity parameter Z max , and the minimum time delay t min d .In Figure 1, it was illustrated that smaller values of Z max , κ and t min d corresponded to larger stochastic gravitational-wave background signals.
A detection of a stochastic background by definition requires a reasonably large stochastic signal, pushing our posteriors on all three parameters to smaller values in Figure 5.
In addition to the posteriors, we show the merger rate dR dm1 (m 1 = 20M ⊙ ) constructed from the posterior samples in Figure 4 (left), as well as the posterior samples on Ω(f ) in Figure 4 (right).The posterior samples of both quantities show agreement with the injected quantities, which are denoted by the red dash-dotted line in both figures.
Case 2: A non-detection -We now repeat the above analysis, but with a gravitational-wave background mea- surement at A+ sensitivity that falls below the expected sensitivity of a gravitational-wave background search, as displayed in Figure 6 (right).This is illustrated by the fact that the injected gravitational-wave background falls below the A+ PI curve.Similarly to the detectable case, the simulated gravitational-wave background measurement is constructed from a posterior sample of the run with current data to ensure consistency with observed GWTC-3 binary black hole events.In Figure 7, we show that the improved A+ sensitivity, even in the absence of a gravitational-wave background detection, allows for better constraints on the time-delay parameters.In particular, we note that the slope κ of the time-delay distribution and the metallicity threshold Z max are now required to lie along a preferred twodimensional contour.The same is true of κ and the minimum time delay t min d .The result is that very large and negative values of κ, and small values of Z max and t min d are excluded, since these values would all have yielded a detectable gravitational-wave background.Larger values of the metallicity parameter Z max and the minimum time delay t min d are favored to be consistent with a nondetection of a gravitational-wave background.Furthermore, we note that dR dm1 (m 1 = 20M ⊙ ) is poorly constrained in the case a gravitational-wave background detection cannot be claimed, as illustrated in the left-hand side of Figure 6.In particular, we point out the similarity of this posterior with the result with current data presented in the left-hand side of Figure 3, but note the somewhat lower 95% upper limit due to the improved A+ sensitivity.Furthermore, the Ω(f ) posterior in Figure 6 (right) is consistent with the results using current data shown in the right-hand side of Figure 3, although with a slightly lower 95% upper limit, again as a result of the improved A+ sensitivity.
It is particularly interesting to stress that distinct regions of parameter space can be constrained depending on whether the gravitational-wave background is detected or not.The non-detection of a gravitational-wave background at A+ sensitivity would allow us to heavily constrain the slope of the time-delay distribution κ, as other values of that parameter would have given rise to a detectable gravitational-wave background, which was not observed.However, the Z max -t min d parameter space remains mostly unconstrained, although slightly favoring larger values of both parameters.In contrast, the detection of a gravitational-wave background instead allows us to place stronger bounds on Z max and t min d , but yields less stringent constraints on the κ parameter, favoring larger negative values.
To illustrate why distinct regions of parameter space can be constrained with the detection or non-detection of a gravitational-wave background, we consider the signal-to-noise ratio given by where T denotes the observation time, γ 12 the overlap reduction function, and P i the noise power spectral density of detector i = {1, 2} (Romano & Cornish 2017).In Figure 8, we show contours of the stochastic signalto-noise ratio in the κ − Z max parameter space, using the Advanced LIGO A+ sensitivity.When computing the expected gravitational-wave background spectrum Ω(f ), the minimum time delay was set to t min d = 0.05 Gyrs and the values of other hyperparameters were set to the median values of the posteriors in Figure 2. Additionally, we show 1D histograms of regions of parameter space where the expected signal-to-noise ratio is smaller (larger) than 3 in blue (green).In particular, we want for the joint analysis using both direct detections and stochastic background constraints at Advanced LIGO A+ sensitivity, in the presence of an undetectable gravitational-wave background.The dashed black lines represent the 1D histograms for the joint analysis on data from the first three observing runs from Figure 2, for reference.We note that large negative values for the slope of the time-delay distribution κ are disfavored, and there is an enhanced support for larger values of the metallicity parameter Z max and the minimum time delay t min d .In addition, we point out the complementarity of these constraints with the case of a detectable gravitational-wave background at A+ sensitivity, as reported in Figure 5.
to highlight that the regions of parameter space where the signal-to-noise ratio is larger (smaller) than 3 corresponds to the region of parameter space where most posterior samples lie for a detectable (undetectable) gravitational-wave background in Figure 4 (Figure 6).The detection (non-detection) effectively forces the posterior samples to lie within a region of parameter space with a corresponding signal-to-noise ratio that is above (below) the detection threshold of ρ = 3.
Therefore, the two cases considered in this work, both a gravitational-wave background detection and non-detection at the Advanced LIGO A+ sensitivity, provide complementary information about the time-delay distribution and metallicity-specific star formation rate.

CONCLUSION
The target of this work is to shed light on the evolution of binary black holes with gravitational-wave data.We currently have two sources of information -the direct binary black hole detections comprised in GWTC-3, and the upper limits on the stochastic gravitationalwave background.Each of these observational inputs In constructing the stochastic signal-to-noise ratio contours, the minimum time delay was set to t min d = 0.05 Gyrs and all other hyperparameters were set to the median value of the posteriors in Figure 2. The regions of parameter space where the signal-to-noise ratio lies above (below) 3 corresponds to the region of parameter space where most posterior samples lie for a detectable (undetectable) gravitational-wave background in Figure 4 (Figure 6).
provides us information about the binary merger rate at different redshifts.Our goal is to synthesize both together to constrain the parameters that govern the metallicity-specific formation rate of binary black hole progenitors and their subsequent evolutionary time delays We consider LIGO-Virgo data from the first three observing runs, including both individual binary black hole mergers from the GWTC-3 catalog, as well as current upper limits on the gravitational-wave background.Current individual binary black hole mergers allow to constrain the slope of the time-delay distribution to negative values, but leave the other parameters of interest unconstrained.The addition of information from the gravitational-wave background in our joint analysis using individual binary black hole mergers and the upper limits on a gravitational-wave background does not provide more stringent constraints at current sensitivity.
Nevertheless, we consider the future Advanced LIGO A+ sensitivity and consider both the case of a gravitational-wave background detection and nondetection at that sensitivity.We show that both cases offer unique and complementary constraints on the parameter space of interest.Indeed, the non-detection of a gravitational-wave background at Advanced LIGO A+ sensitivity results in tight constraints on the slope of the time-delay distribution.On the other hand, the detection of a gravitational-wave background would allow to constrain larger parameter space regions of the metallicity threshold and minimum time-delay parameter.
Although considering current gravitational-wave background upper limits does not improve the constraints on the metallicity-specific star formation rate and time-delay parameters at current sensitivity, this work shows that the additional information contained in a gravitational-wave background measurement will be essential to learn more about the environment in which compact binaries formed and their evolution.yield distinct and complementary information about the environment in which binaries formed and their evolution throughout time.
In addition, our choice of model for the cumulative metallicity of star formation, as given by Eq. ( 2), represents a similar systematic source of uncertainty.Future work will explore how variations in metallicity models similarly impact our conclusions.for the joint analysis using both direct detections and stochastic background constraints at Advanced LIGO A+ sensitivity, in the presence of a detectable gravitational-wave background.The dashed black lines represent the 1D histograms for the joint analysis on data from the first three observing runs, for reference.We note that there is more support for larger negative values of the slope of the time-delay distribution κ, and smaller values of the metallicity parameter Z max and the minimum time delay t min d are favored.In addition, we point out the complementarity of these constraints with the case of an undetectable gravitational-wave background at A+ sensitivity, as reported in Figure 11.

B. ADDITIONAL DISTRIBUTIONS AND PRIORS
In addition to the mass and redshift distribution described in the main text, we provide the assumed distributions for the spin magnitudes and spin-orbit tilt angles.These are modeled as truncated Gaussian distributions (Callister where the mean is assumed to be 1, and the variance σ 2 u is inferred from the data.Furthermore, we provide a list of priors used for the parameter estimation performed in this work in Table 1.  1. Prior choice for the hyperparameters describing the mass distributions, the time-delay distribution and the spin distribution.Note that, for the minimum time-delay parameter t min d , we do not extend the prior to lower values to not be in conflict with the minimum lifetime of massive stars (e.g., Kippenhahn et al. 2013).

C. DATA
In this appendix, we provide more detailed information regarding the exact inputs to our analysis.In our analyses, we include the binary black holes of the GWTC-3 catalog (Abbott et al. 2023a) detected with a false-alarm rate below 1 yr −1 .In the GWTC-3 catalog, there are two events, GW190814 (Abbott et al. 2020) and GW190917, that are possibly binary black holes but which are known to be outliers with respect to the bulk binary population (Abbott et al. 2023a).We do not consider these two events, leaving 69 binary black holes to be included in our analysis.Publiclyavailable parameter estimation samples are used, which are provided by the Gravitational-Wave Open Science Center3 (Vallisneri et al. 2015;Abbott et al. 2021dAbbott et al. , 2023b) ) and/or Zenodo.We use the "Overall posterior" samples4 for GWTC-1 events (Abbott et al. 2019), use the "PrecessingSpinIMRHM" samples for events from GWTC-2 (Abbott et al. 2021c)5 , and use use the "C01:Mixed" samples6 for new events in GWTC-3 (Abbott et al. 2023a).These contain the combination of parameter estimation samples from several waveform families, each of which include the physical effects of spin precession.We note that the GWTC-2 and GWTC-3 samples additionally include higher-order mode content, whereas GWTC-1 samples, in contrast, do not.
In our analysis, injected signals are used to characterize selection effects, as given by Eq. ( 18).We use the injection set7 detailed in Abbott et al. (2023a), labeling injections as "found" when they are recovered with false-alarm rates below 1 yr −1 in at least one search pipeline.However, the injections for the O1 and O2 observing runs do not have associated false-alarm rates, only network signal-to-noise ratios ρ.We therefore consider them "found" if ρ ≥ 10.
In addition, for the joint analysis using both direct detections and stochastic background constraints at current O3 sensitivity, we use the results of the search for an isotropic gravitational-wave background following the LIGO-Virgo O3 observing run8 (Abbott et al. 2021a).Note that we only consider the LIGO Hanford -LIGO Livingston baseline, since the LIGO-Virgo baselines have negligible sensitivity to CBC signals.

D. COMPUTATION OF Ω(f )
The computation of the expected Ω(f ), given a set of hyperparameters Λ, can be written as in Eq. ( 6) in the main text.We recall that this depends on the average energy dEs dfs radiated by each binary, given by Eq. ( 7), where the integration is performed over the masses m 1 and m 2 .Therefore, including the integration over redshift, the computation of the expected Ω(f ) requires the evaluation of a three-dimensional integral.Numerically, the integral can be evaluated on a grid, resulting in a fairly inefficient and slow implementation of the computation.Alternatively, one can rely on a Monte Carlo integration technique to evaluate the integral, where one relies on averaging many individual draws of the population.
In our implementation of this approach, N samples are drawn from a uniform distribution at the start of the computation for each of the parameters z, m 1 , and m 2 .We denote the i-th parameter draw as ϕ i = {z i , m i 1 , m i 2 }.We proceed by computing the energy spectra dE s /df s for each of these parameter draws and label these dE i s /df s .However, we want to compute the Ω(f ) for the population corresponding to the distributions described in the main text, with hyperparameters Λ.We therefore rely on reweighting dE i s /df s given some hyperparameters Λ. Schematically, the computation takes the form where the average is taken over the number of sample draws N , and the weights w i are defined as p(m i 2 ) p draw (m i 2 ) .

(D5)
The p draw distributions represent the uniform distributions used to make the initial draws of the parameters ϕ i , whereas the numerators denote the distributions we are reweighting to, characterizing the population, corresponding to a set of hyperparameters Λ, for which we want to compute the Ω(f ) spectrum.
In order to show the convergence of the Monte Carlo average, we show the Ω(f ) spectrum computed by performing the integral explicitly and compare it to our implementation for several values of the number of initial draws.The results are displayed in Figure 12.Note that although the results are only shown for one set of hyperparameters, i.e., one specific distribution, and therefore, one Ω(f ) spectrum in the plot, tests were performed on various other distributions to verify convergence.In general, a value N = 20000 approximates the full integral computation well in the frequency range of interest for gravitational-wave background searches with the LIGO-Virgo-KAGRA detectors, with relative errors of O(10 −1 − 10 −2 ).Therefore, throughout our work, we use N = 20000 when computing the Ω(f ) spectrum with the Monte Carlo averaging procedure.

Figure 1 .
Figure 1.Left: Examples of the binary black hole merger rate R(z) for different values of the slope of the time-delay distribution, κ, and different values of the maximum metallicity below which black holes are formed, Z max .All curves are arbitrarily normalized to R(z = 0.2) = 1.The value of the minimum time-delay parameter, t min d , is set to 0.05 Gyrs.Steeper time-delay distributions result in a larger merger rate at large redshifts, whereas larger metallicity thresholds result in a lower merger rate at large redshifts.Right: Example of Ω(f ) spectrum for different values of the same parameters κ and Z max .The value of all other parameters that enter in the mass distribution as given by Eqs.(14) -(16) is fixed to arbitrary values: m low = 10M⊙, m high = 80M⊙, α = −2, µm = 20M⊙, σm = 3M⊙, fp = 0.01, βq = 2, R ref = 1 M −1 ⊙ Gpc −3 yr −1 , δm low = 1M⊙, δm high = 3M⊙, and t min Figure 1.Left: Examples of the binary black hole merger rate R(z) for different values of the slope of the time-delay distribution, κ, and different values of the maximum metallicity below which black holes are formed, Z max .All curves are arbitrarily normalized to R(z = 0.2) = 1.The value of the minimum time-delay parameter, t min d , is set to 0.05 Gyrs.Steeper time-delay distributions result in a larger merger rate at large redshifts, whereas larger metallicity thresholds result in a lower merger rate at large redshifts.Right: Example of Ω(f ) spectrum for different values of the same parameters κ and Z max .The value of all other parameters that enter in the mass distribution as given by Eqs.(14) -(16) is fixed to arbitrary values:m low = 10M⊙, m high = 80M⊙, α = −2, µm = 20M⊙, σm = 3M⊙, fp = 0.01, βq = 2, R ref = 1 M −1 ⊙ Gpc −3 yr −1 , δm low = 1M⊙, δm high = 3M⊙,and t min d = 0.05 Gyrs.For more details on the time-delay parameter t min d and R ref , we refer the reader to Eqs.(3) and (12), respectively.The behavior of the Ω(f ) spectrum as the time-delay parameters are varied can directly be related to the behavior of the merger rate on the left-hand side, as a larger integrated merger rate results in a larger net Ω(f ) spectrum.

Figure 2 .
Figure 2. Posteriors on the parameters governing binary black hole birth and evolution using both direct detections and stochastic background constraints on the first three LIGO-Virgo observing runs.For clarity, the only posteriors that are shown are for the reference merger rate amplitude R ref , the slope of the time-delay distribution κ, the metallicity threshold Z max , and the minimum time delay t min d , although we simultaneously inferred the binary black hole mass distribution as in Sec. 3. The reference merger rate amplitude R ref is well constrained due to the currently observed individual binary black hole merger events happening at low redshifts.Although the metallicity threshold Z max and the minimum time-delay parameter t min d cannot be constrained with data from the first three observing runs, the slope of the time-delay distribution κ is constrained to negative values.Positive values of this parameter would give rise to a different merger rate that decreases with redshift, contradicting current observations.
Figure3.Left: Posterior samples on the merger rate dR dm 1 (m1 = 20M⊙) for the joint analysis using both direct detections and stochastic background constraints.For reference, the 5% and 95% credible limits are shown in full black lines for the individual binary black hole run, and the dashed black line shows the median value for the binary black hole only run.The results of the individual binary black hole run and the joint analysis using both direct detections and stochastic background constraints are comparable, illustrating that the gravitational-wave background upper limits are not informative at current detector sensitivity.Furthermore, we note that the merger rate is well-constrained at low redshifts because of the individual binary black hole mergers detected at those redshifts.Right: Posterior on Ω(f ) using the posterior samples from the joint analysis using both direct detections and stochastic background constraints on the first three LIGO-Virgo-KAGRA observing runs.The median is denoted by the dashed black line, and the 95% credible upper limit by the full black line.The 2σ O3 PI curve fromAbbott et al. (2021a) is shown (dotted black) as an indication of the gravitational-wave background detection threshold after the first three LIGO-Virgo observing runs.All posterior samples fall below the O3 PI curve, illustrating that we do not expect the gravitational-wave background contribution to be informative at current detector sensitivity.

Figure 4 .
Figure 4. Left: Posterior samples for dR dm 1 (m1 = 20M⊙) with the injected rate corresponding to a detectable gravitational-wave background at A+ sensitivity indicated by the red dash-dotted line.The full black lines denote the 5-95% confidence interval, and the dashed black line denotes the median.Right: Posterior samples for Ω(f ) with the injected detectable background denoted by the red dash-dotted line.The 95% upper limit is denoted by the full black line, and the dashed black line represents the median.The dotted black line represents the 2σ A+ PI curve.

Figure 5 .
Figure 5. Posterior on the slope of the time-delay distribution κ, the metallicity threshold Z max , and the minimum time-delay parameter t min d

Figure 6 .
Figure6.As in Figure4, but now for the case of an undetectable gravitational-wave background.

Figure 7 .
Figure 7. Posterior on the slope of the time-delay distribution κ, the metallicity threshold Z max , and the minimum time-delay parameter t min d

Figure 8 .
Figure8.Contours of the expected stochastic signal-tonoise ratio using the Advanced A+ sensitivity for the slope of the time-delay distribution κ and the maximum metallicity below which black holes are formed Z max .The 1D histograms denote the regions of parameter space where the signal-to-noise ratio is smaller (larger) than 3 in blue (green).In constructing the stochastic signal-to-noise ratio contours, the minimum time delay was set to t min

Figure 10 .
Figure 10.Posterior on the slope of the time-delay distribution κ, the metallicity threshold Z max , and the minimum time-delay parameter t min d