Shouts and Murmurs: Combining Individual Gravitational-Wave Sources with the Stochastic Background to Measure the History of Binary Black Hole Mergers
Thomas Callister Maya Fishbach Daniel Holz Will Farr
astro-ph.HE gr-qc



Tom Callister
tcallister@flatironinstitute.org





Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA


Enrico Fermi Institute, Department of Physics, Department of Astronomy and Astrophysics, and Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA







One of the goals of gravitational-wave astronomy is to quantify the evolution of the compact binary merger rate with redshift.
The redshift distribution of black hole mergers would offer considerable information about their evolutionary history, including their progenitor formation rate, the dependence of black hole formation on stellar metallicity, and the time delay distribution between formation and merger.
Efforts to measure the binary redshift distribution are currently limited, however, by the detection range of existing instruments, which can individually resolve compact binary merger events only out to z≲1.
We present a novel strategy with which to measure the redshift distribution of binary black hole mergers well beyond the detection range of current instruments.
By synthesizing direct detections of individually resolved mergers with indirect searches for the stochastic gravitational-wave background due to unresolved distant sources, we can glean information about the peak redshift, z_p, at which the binary black hole merger rate attains its maximum, even when this redshift is beyond the detection horizon.
Using data from Advanced LIGO and Virgo's first and second observing runs, we employ this strategy to place joint constraints on z_p and the slope α with which the binary merger rate increases at low redshifts, ruling out merger rates that grow faster than α≳ 7 and peak beyond z_p ≳ 1.5.
Looking ahead, we project that approximately one year of observation with design-sensitivity Advanced LIGO will further break remaining degeneracies, enabling a direct measurement of the peak redshift of the binary black hole merger history.







§ INTRODUCTION



The Advanced LIGO <cit.> and Advanced Virgo <cit.> gravitational-wave experiments are rapidly transitioning between low- and high-statistics regimes.
With the LIGO-Virgo detections of eleven compact binary mergers during the past O1 and O2 observing runs <cit.> and tens more anticipated in the present O3 run <cit.>, we can now begin to understand the ensemble properties of compact binaries, including the distributions of their component masses and spins  <cit.>.
Beyond the distributions of these intrinsic binary parameters, we might also seek to understand the redshift distribution of binary black hole mergers – how the merger rate evolves as we look back to earlier times in the Universe's history.
If measured, the redshift distribution of compact binary mergers would offer substantial insight into the birth and evolution of compact binaries, encoding such properties as the time delay distribution between black hole formation and merger <cit.>, the dependence of black hole production on stellar metallicity <cit.>, and perhaps even the relative contributions from competing binary formation channels; e.g. field binaries, hierarchical triples, dynamical capture, or primordial black holes <cit.>.

Study of the binary black hole redshift distribution, however, is made difficult by the limited range of existing gravitational-wave detectors.
Figure <ref>, for example, shows a typical model for the source-frame rate of binary black hole mergers as a function of redshift.
To obtain this figure, we assume progenitor formation following the star formation rate of <cit.> weighted by the fraction of stellar formation occurring at metallicities Z≤ 0.3 Z_⊙ <cit.>.
We further adopt a p(t_d) ∝ t_d^-1 probability distribution for the time delay t_d between binary formation and merger, with 50 Myr≤ t_d ≤ 13.5 Gyr.
Within this simple model, the binary black hole merger rate peaks at z∼2, while more sophisticated models generally predict merger rates peaking between redshifts z∼2 to 4, depending on the specific formation channel presumed <cit.>.

In contrast, design-sensitivity Advanced LIGO is expected to successfully detect optimally-oriented 30+30 M_⊙ binary black holes only out to redshifts z≲1.2 <cit.>.
Current efforts to study the redshift distribution of compact binary mergers therefore attempt only to measure the leading-order, low-redshift evolution of the binary merger rate <cit.>; observation of the peak and subsequent turnover of the black hole redshift distribution is a challenge left to future third-generation detectors <cit.>.

In this paper we demonstrate that present-day gravitational-wave observatories can provide meaningful measurements of the high-redshift evolution of the compact binary merger rate.
We achieve these measurements by synthesizing the direct detections of compact binaries in the local Universe with an additional piece of information: the astrophysical stochastic gravitational-wave background <cit.>.
Arising from the superposition of all distant individually-undetectable compact binaries, the stochastic gravitational-wave background manifests as excess correlated noise shared amongst a network of gravitational-wave detectors <cit.>.
The strength of the present-day gravitational-wave background is determined by the cumulative merger history of binary black holes, integrated across all redshifts <cit.>.
The observation of (or even upper limits on) the gravitational-wave background can therefore be leveraged to place powerful constraints on the redshift distribution of binary mergers, complementary to those constraints gleaned from the direct detection of binaries in the local Universe (Sec. <ref>).

We apply our approach to existing data, finding that the synthesis of binary black hole detections <cit.> and gravitational-wave background constraints <cit.> from Advanced LIGO and Advanced Virgo's first two observing runs already yields non-trivial constraints on the peak of the binary black hole redshift distribution (Sec. <ref>).
With additional data gathered from future observing runs, our method may enable a measurement of this peak redshift within the next five years (Sec. <ref>).




§ HIGH REDSHIFT CONSTRAINTS FROM THE GRAVITATIONAL-WAVE BACKGROUND






In their O1 and O2 observing runs, Advanced LIGO and Virgo confidently detected ten binary black hole mergers, the most distant of which (GW170729) may have occurred at z≈0.5 <cit.>.
Together, these ten events have recently allowed for the first exploration of the binary black hole merger rate's evolution with redshift.
Adopting a model
	
    ℛ(z) = ℛ_0 (1+z)^α,

for the source-frame merger rate per comoving volume <cit.>, <cit.> find α = 6.5^+9.1_-9.3 at 90% credibility.
Thus, in the local universe, the binary black hole merger rate (probably) increases with redshift.

If the binary black holes observed with LIGO and Virgo are born from stellar progenitors, then the black hole merger rate cannot continue to increase out to arbitrarily high redshifts.
Instead, it must reach a maximum at some peak redshift, z_p, and then decay to zero as star formation ceases in the very early Universe.
Generically, we can describe this complete merger history with a phenomenological model of the form <cit.>:
	
    ℛ(z) = 𝒞(α,β,z_p) ℛ_0  (1+z)^α/1+(1+z/1+z_p)^α+β ,

allowing a source-frame merger rate that initially evolves as ℛ(z) ∝ (1+z)^α, reaches a maximum near z_p, and subsequently falls as ℛ(z)∝ (1+z)^-β.
The example binary black hole merger rate shown in Fig. <ref>, for example, is well-fit by this phenomenological model using α = 1.9, β = 3.4, and z_p = 2.4, shown via a dashed grey curve.
The normalization constant 𝒞(α,β,z_p) = 1+(1+z_p)^-α-β ensures that ℛ(0) = ℛ_0.

At present the direct detection of binary black holes with Advanced LIGO and Virgo can offer no meaningful constraints on z_p or β.
In O2, the range within which Advanced LIGO could detect a typical 30+30 M_⊙ binary black hole <cit.> was z≲0.5; in the future O5 observing run this range may be pushed to z≲1.2 <cit.>.
Meanwhile, if the black hole merger rate roughly follows the star formation rate, it should peak at z_p≳2, well beyond our ability to probe with direct detections.

We have another piece of information at our disposal, however.
Although individually undetectable, the superposition of all distant binary black holes gives rise to a stochastic gravitational-wave background, detectable in the form of excess cross-power between widely-separated detectors <cit.>.
The stochastic gravitational-wave background is conventionally described by a dimensionless energy-density spectrum <cit.>
	
    Ω(f) = 1/ρ_cdρ_gw/dln f,

where dρ_gw/dln f is the present-day energy density in gravitational-waves per logarithmic frequency interval and ρ_c = 3 H_0^2 c^2/8π G is the Universe's critical energy density.
Here, c is the speed of light, G is Newton's constant, and H_0 is Hubble's constant; we adopt H_0 = 70 km s^-1 Mpc^-1.

The energy density arising from the population of binary black hole mergers is given by <cit.>:
	
    Ω(f) = f/ρ_c∫_0^z_max dz ℛ(z) ⟨dE_s/df_s|_f(1+z)⟩/(1+z) H(z).

Here, ⟨ dE_s/df_s ⟩ is the source-frame energy spectrum radiated by a single binary <cit.>, averaged over the binary black hole population.
If the intrinsic parameters of individual binary black holes (e.g. their masses and spins) are denoted by ϕ and have distribution p(ϕ), then
	
    ⟨dE_s/df_s⟩ = ∫ dϕ p(ϕ) dE_s/df_s(ϕ).

Note that in Eq. (<ref>) we evaluate ⟨ dE_s/df_s ⟩ at the source-frame frequency f(1+z).
Meanwhile, H(z) = H_0 √(Ω_m(1+z)^3 + Ω_Λ) is the Hubble parameter at redshift z (neglecting radiation density).
We take the energy densities of matter and dark energy to be Ω_m = 0.3 and Ω_Λ = 0.7, respectively.
Finally, the integral in Eq. (<ref>) is taken up to a cutoff redshift z_max; we fix z_max = 10, beyond which we expect virtually no star formation and hence no black hole mergers (assuming stellar progenitors).
Alternatively, allowing z_max itself to vary as another free parameter may help to provide constraints on binary black holes of non-stellar origin, like the mergers of primordial black holes <cit.>.

The energy density, Ω(f), measured by stochastic searches is, in essence, a weighted integral over the binary black hole merger history ℛ(z), sensitive to the total number of past  mergers.
Thus, if the local rate ℛ_0 is independently fixed by direct detections, then knowledge of Ω(f) provides strong bounds on the possible values of α, β, and z_p.
This is true even given a non-detection of the gravitational-wave background.
To illustrate this, we can consider how the signal-to-noise ratio (S/N) of the gravitational-wave background varies with α and z_p.

Given a model Ω_M(f) for the true energy-density spectrum, the signal-to-noise ratio of the gravitational-wave background is <cit.>
	
    S/N = (Ĉ | γ Ω_M)/√((γ Ω_M | γ Ω_M)).

Here,
	
    Ĉ(f) = 1/T20π^2/3 H_0^2 f^3 s̃_1(f) s̃^*_2(f)

is the cross-correlation statistic between the strains s̃_1(f) and s̃_2(f) measured by two gravitational-wave detectors <cit.>, and we have defined an inner product
	
    (A|B) = 2 T (3H_0^2/10π^2)^2 ∫_0^∞ df Ã(f) B̃^*(f)/f^6 P_1(f) P_2(f),

where P_i(f) is the one-sided noise power spectral density of detector i and T is the total observation time.
In the presence of a gravitational-wave background, the expectation value of Ĉ(f) is
	
    ⟨Ĉ(f) ⟩ = γ(f) Ω(f)

and its variance is ⟨Ĉ(f) Ĉ(f')⟩ = δ(f-f')σ^2(f), with
	
    σ^2(f) = 1/T(10π^2/3 H_0^2)^2 f^6 P_1(f) P_2(f).

The factor γ(f) in Eqs. (<ref>) and (<ref>), known as the overlap reduction function, quantifies the geometrical sensitivity of a given detector pair to an isotropic gravitational-wave background <cit.>.
The optimal S/N occurs when we choose a model Ω_M(f) = Ω(f) matching the true energy density; the expected S/N in this case is <cit.>
	
    ⟨S/N⟩_opt  = √((γ Ω | γ Ω)).




In Fig. <ref> we plot the optimal S/N with which the gravitational-wave background would have appeared in O1 and O2 as a function of possible values for α and z_p.
In this example we fix ℛ_0 = 30 Gpc^-3 yr^-1 and β = 3, and assume a population of equal mass binaries with chirp mass ℳ_c = 30 M_⊙.
If α≲5, virtually no stochastic signal is expected, consistent with the non-detection of the gravitational-wave background in O1 <cit.> and O2 <cit.>.
However, the expected S/N rises sharply towards the upper-right corner of Fig. <ref>.
In particular, if α≳ 5 and z_p≳1, we should have seen an extraordinarily loud stochastic gravitational-wave signal.
The fact that no such background was detected means that we can already reject this portion of parameter space, ruling out binary black hole backgrounds rising faster than α∼5 and peaking beyond z_p∼1.
We note, though, that these exact limits depend strongly on the assumed local merger rate ℛ_0 and black hole mass distribution (and to a lesser extent on β), and so the results in Fig. <ref> should be taken as an example only.
In Sect. <ref> below, we will instead seek to simultaneously measure these different properties, leveraging both the observational limits on the stochastic gravitational-wave background and the current catalog of direct binary black hole detections.

So far, our argument has implicitly assumed that the distribution of binary black hole parameters is independent of redshift, such that the average energy radiated by a given binary [Eq. (<ref>)] does not vary with z.
This is not necessarily the case.
It is possible, for instance, that black holes born at high redshifts are preferentially more massive, due to the increased stellar masses predicted to occur at low metallicities <cit.>, although more recent work suggests that the mass distribution of merging binaries may be approximately constant <cit.>.
By neglecting the possibility increased masses at higher redshifts, the constraints we obtain on α and z_p are conservative.
Given a fixed observational limit on Ω(f), any presumed increase in the average radiated energy ⟨ dE_s/df_s ⟩ must be balanced be a decrease in the merger rate ℛ(z) at high redshifts, yielding stricter limits than those shown in Fig. <ref>.
Nevertheless, one could incorporate effects like metallicity-dependent masses in this analysis by amending Eq. (<ref>) to additionally include integration over distributions of formation redshifts or progenitor metallicities <cit.>.




§ PEAK REDSHIFT CONSTRAINTS FROM O1 AND O2




The best constraints on ℛ(z) will come from neither the direct detection of binary black holes nor the gravitational-wave background searches considered separately, but instead from a joint analysis that self-consistently synthesizes both sources of information.
In this paper we perform the first such joint analysis, synthesizing stochastic data and direct black hole observations to hierarchically measure the redshift distribution of binary black hole mergers.
We take as inputs the integrated cross-correlation spectrum Ĉ(f) measured between the LIGO Hanford and Livingston detectors <cit.> during O1 and O2, as well as parameter estimation results for each of the ten binary black hole mergers comprising the LIGO and Virgo GWTC-1 catalog <cit.>.

In order to robustly constrain ℛ(z), it will also be important to simultaneously fit for the mass distribution of binary black holes.
In Fig. <ref> above, for instance, the exact exclusion region depends on our presumed black hole mass distribution: heavier or lighter black holes would increase or decrease the expected energy density Ω(f), leading us to draw different conclusions about ℛ(z) in the case of a stochastic non-detection.
Strong degeneracies also exist between the inferred mass and redshift distributions of directly-detected black hole mergers <cit.>.
A dearth of detections at large redshifts, for example,  simply implies a low merger rate for high-mass binaries, since low-mass binaries go undetected at large distances.
This can be explained either by a low overall rate at high redshifts, or by a mass distribution that prefers low-mass binaries.

Consider a population of binary black hole mergers, with a local merger rate per unit comoving volume ℛ_0 and whose mass and redshift distributions are characterized by parameters Λ.
The likelihood of obtaining data {d_i}_i=1^N_obs from N_obs direct detections, as well as a stochastic cross-correlation spectrum Ĉ(f), is
	
    p(Ĉ,{d_i} | Λ,ℛ_0) =  p_bbh({d_i} | Λ,ℛ_0) p_stoch(Ĉ | Λ,ℛ_0),

which has been factored into a direct-detection and a stochastic term.

The likelihood, p_bbh({d_i} | Λ,ℛ_0), of our direct binary black hole detections is given by <cit.>
	
    p_bbh({d_i}|Λ,ℛ_0) 
       ∝[N(Λ,ℛ_0) ξ(Λ)]^N_obs e^-N(Λ,ℛ_0)ξ(Λ)
       ×∏_i=1^N_obs∫ p(d_i | ϕ) p(ϕ |Λ) dϕ/ξ(Λ).

Here, p(d_i | ϕ ) is the likelihood for event i given its component masses m_1 and m_2 and redshift z, together abbreviated as ϕ = {m_1,m_2,z}.
Meanwhile, p(ϕ | Λ) is the ensemble distribution of these source parameters.
The quantity N(Λ,ℛ_0) is the total number of binary black hole mergers (both observed and unobserved) expected to occur during our observation time; see Eq. (<ref>) below.
Observational selection effects are captured by the factor ξ(Λ), the fraction of all binary black holes that we expect to successfully detect.
If P_det(ϕ) is the probability of successfully detecting an event with parameters ϕ, then
	
    ξ(Λ) = ∫ P_det(ϕ) p(ϕ|Λ) dϕ.

In our analysis we precompute P_det(ϕ) over a grid of masses and redshifts, using the semi-analytic prescription of <cit.>, and requiring detections to have a matched filter signal-to-noise ratio of ρ>8 in a single detector.

In practice, we do not have direct access to the likelihoods, p(d_i|ϕ), needed to compute Eq. (<ref>).
Instead, we have discrete samples {ϕ_i} drawn from each event's posterior distribution p(ϕ | d_i), obtained via parameter estimation with Monte Carlo integration or nested sampling <cit.>.
Parameter estimation itself is performed while assuming some default prior, p_pe(ϕ), that is generally not equal to the population prior p(ϕ|Λ) appearing in Eq. (<ref>).
To evaluate Eq. (<ref>), we must therefore replace the integral with an average over discrete samples, weighting each sample with p_pe^-1(ϕ) to undo the influence of the prior used in parameter estimation:
	
    p_bbh({d_i}|Λ,ℛ_0) 
       ∝[N(Λ,ℛ_0) ξ(Λ)]^N_obs e^-N(Λ,ℛ_0) ξ(Λ)
       ×∏_i=1^N_obs1/ξ(Λ)⟨p(ϕ_i | Λ)/p_pe(ϕ_i)⟩_samples.


The stochastic cross-correlation spectrum Ĉ(f), meanwhile, is generally obtained through the weighted combination of a large number of measurements performed over short 𝒪(100 s) time segments <cit.>, and so the likelihood p_stoch(Ĉ |Λ,ℛ_0) is well-approximated as a Gaussian <cit.>:
	
    p_stoch(Ĉ |Λ,ℛ_0) 
       ∝exp[- 1/2(Ĉ - γ Ω_M(Λ,ℛ_0) | Ĉ - γ Ω_M(Λ,ℛ_0))],

where Ω_M(Λ,ℛ_0;f) is our model energy-density spectrum and we have used the inner product defined in Eq. (<ref>).

We model the intrinsic redshift distribution of binary black hole mergers as
	
    p(z|α,β,z_p) ∝1/1+zℛ(α,β,z_p;z) dV_c/dz,

where ℛ(α,β,z_p;z) is given in Eq. (<ref>) and dV_c/dz is the comoving volume per unit redshift;
note that Eq. (<ref>), once normalized, is independent of the local merger rate ℛ_0.
The leading factor of (1+z)^-1 transforms between source-frame and detector-frame times.
Correspondingly, the total number N of mergers expected to occur during our observation time T is
	
    N(α,β,z_p,ℛ_0) = T ∫_0^z_maxdz  1/1+zℛ(α,β,z_p,ℛ_0;z) dV_c/dz.


Following <cit.> and <cit.>, we model the ensemble distribution of primary black hole masses as a power law
	
    p(m_1 | κ, M_min, M_max) ∝
    	m_1^-κ    (M_min≤ m_1 ≤ M_max) 
    
    	0     (else)

and assume a flat distribution
	
    p(m_2 | m_1, M_min) =
    	1/m_1 - M_min    (M_min≤ m_2≤ m_1) 
    
    	0     (else)

of secondary masses.







In our analysis we fix M_min = 5 M_⊙, while hierarchically inferring the parameters {α,β,z_p,ℛ_0,κ,M_max} of the binary black hole redshift and mass distributions.
We adopt the priors listed in Table <ref>, and perform inference using  <cit.>.
For every iteration of our sampler, we evaluate the direct-detection likelihood in Eq. (<ref>), using Eq. (<ref>) to convert the proposed event rate density ℛ_0 to a total number of mergers N.
We then compute a model stochastic energy-density spectrum, integrating over the proposed mass and redshift distributions [in Eqs. (<ref>) and Eq. (<ref>), respectively] of the binary black hole population, thereby evaluating the stochastic contribution [Eq. (<ref>)] to the overall likelihood.

Figure <ref> shows our resulting posterior on the rate evolution of binary black hole mergers, using the direct GWTC-1 detections alone (left) and combining direct detections with existing stochastic search results (right).
Each trace in these figures represents a draw from our {α,β,z_p,ℛ_0} posterior. The left panel of Fig. <ref> is directly comparable to Fig. 6 of <cit.>.
Figures <ref> and <ref> show the corresponding posteriors on these parameters, marginalized over κ and M_max.
Full parameter estimation results are listed in Table <ref>.

Direct detections alone allow a measurement of the local merger rate to ℛ_0 = ^+_- Gpc^-3 yr^-1 at 95% credibility (the most precise measurement actually occurs at the “waist” seen at z∼0.1). This is consistent with the results of <cit.>.
Direct observations also allow us to roughly constrain α, with a moderate preference for α∼ 5 shown in Fig. <ref>.
Significant uncertainties remain, however.
At 95% credibility, we find α = ^+_-, and, since the α posterior extends all the way to our lower prior bound, we can only robustly constrain α≤.
Direct detections offer no information about β or z_p.
Correspondingly, in Fig. <ref> we have virtually no constraints on the merger rate beyond z∼1.
At z=1.5, for example, the local merger rate could plausibly lie anywhere between 10^-4–10^8 Gpc^-3 yr^-1, a range spanning twelve orders of magnitude.



In contrast, the inclusion of O1 and O2 stochastic search data provides a hard upper bound on the high-redshift merger rate; our non-detection of the gravitational-wave background effectively excludes rate densities that rise above ℛ(z)≳10^4 Gpc^-3 yr^-1.
This additional constraint is reflected in Fig. <ref>.
While the inclusion of O1 and O2 stochastic data does not notably alter the one-dimensional marginal posteriors, it does significantly alter our joint posterior on α and z_p.
As argued in Sec. <ref>, the non-detection of a stochastic gravitational-wave signal provides a joint constraint on these two parameters, rejecting a large portion of the α-z_p parameter space.
When this stochastic exclusion region is combined with the constraint on α from direct GWTC-1 detections, we can already see hints of a preferred contour in the α–z_p plane.

Although the primary goal of this analysis is to measure the evolution of the binary black hole merger rate, it additionally provides a self-consistent framework for predicting the energy density Ω(f) of the binary black hole background using both the known population properties of local binary black holes and upper limits from Advanced LIGO and Advanced Virgo stochastic searches <cit.>.
For every posterior sample in Fig. <ref> (including the mass parameters M_max and κ not shown there) we compute the corresponding binary black hole energy density using Eq. (<ref>).
The result, shown in Fig. <ref>, is a prediction for the binary black hole stochastic background that is marginalized over our uncertainty in both the mass distribution and rate evolution of binary black holes, and subject to the measured upper limits from Advanced LIGO.

Within Fig. <ref>, the dashed black curve traces the 95% credible upper limit on Ω(f) at each frequency.
For comparison, the dashed blue curve shows the 2σ “power-law integrated (PI) curve” <cit.> quantifying Advanced LIGO's integrated sensitivity to the gravitational-wave background following O1 and O2; energy-density spectra lying above this curve will generally be observed with S/N≥ 2.
As expected, the 95% credible limit on Ω(f) lies nearly tangent to the PI curve.
The solid black curve, meanwhile, marks the median predicted energy-density.
At 25 Hz, this median prediction gives Ω(25 Hz) = 8.8×10^-10, comparable to the prediction made by <cit.>: Ω(25 Hz) = 5.3×10^-10.
The uncertainty on our predicted energy-density spectrum, though, is considerably larger.
While the <cit.> model includes uncertainty on the local rate density ℛ_0 of binary black hole mergers, it makes stringent assumptions concerning the subsequent evolution of the merger rate with redshift, assumptions that carry considerable systematic uncertainty.
In contrast, Fig. <ref> includes marginalization over all possible redshift distributions, making this systematic uncertainty explicit.




§ ADVANCED LIGO AT DESIGN SENSITIVITY




The continued synthesis of direct detections with stochastic search results will offer increasingly strong information regarding the leading slope, α, and peak, z_p, of the binary black hole merger history.
Additional binary black holes detected in the local Universe will yield ever tighter posteriors on α, while continued time integration by stochastic searches will reject a growing fraction of the joint α–z_p posterior space.
Eventually these two effects will meet, converging to produce a true measurement of both α and z_p.

To illustrate this, here we anticipate the results that will soon be possible with design-sensitivity Advanced LIGO.
We simulate a mock catalog of 500 binary black hole detections, drawn from a population whose mass distribution is characterized by κ=1.2, M_max =45 M_⊙, and M_min=5 M_⊙.
We assume a redshift distribution given by α=3, β=3, z_p=2, and ℛ_0=30 Gpc^-3 yr^-1.
With this choice of local merger rate, we would expect to detect these 500 binary black holes after T∼1.2 years of observation with design-sensitivity Advanced LIGO.




We construct our mock catalogue following <cit.>.
For each event, we draw an “observed” signal-to-noise ratio
	
    ρ_obs∼𝒩(ρ,1)

from a Gaussian distribution about the event's true signal-to-noise ratio ρ, calculated in a detector with a noise power spectral density given by the Advanced LIGO “design sensitivity" curve of <cit.>.
We require our detected events to have ρ_obs>8 in a single detector.
For each detected event, we draw an observed maximum-likelihood chirp mass
	
    logℳ_obs∼𝒩(logℳ ,σ_ℳ8/ρ_obs)

and symmetric mass ratio
	
    η_obs∼𝒩(η ,σ_η8/ρ_obs),

where ℳ and η are the event's true parameters and we adopt characteristic uncertainties σ_ℳ = 0.08 and σ_η = 0.022.
We then draw synthetic likelihood samples about logℳ_obs and η_obs, with variances consistent with the above distributions. This prescription gives realistic uncertainties on the measured component masses and distances of BBH detections, matching the typical uncertainties reported in <cit.>.

We encapsulate a binary's inclination angle and sky location in a single <cit.> projection factor Θ, which quantifies a signal's amplitude reduction due to suboptimal viewing angles and/or sky placement.
If ρ_opt is a binary's optimal signal-to-noise ratio (i.e. face-on and directly overhead), then Θρ_opt is the event's actual signal-to-noise ratio.
For each mock event, we draw a maximum-likelihood projection factor from
	
    Θ_obs∼𝒩(Θ,σ_Θ8/ρ_obs)

where σ_Θ = 0.15, about which we draw likelihood samples {Θ}.

Realistic redshift samples will be strongly correlated with an event's recovered S/N as well as its projection factor Θ.
To capture these correlations, we first draw S/N ratio samples
	
    {ρ}∼𝒩(ρ_obs,1).

Then, noting that ρ is inversely proportional to an event's luminosity distance D_L, we convert {ρ} and {Θ} into luminosity distance samples via
	
    {D_L}/1 Gpc =  ρ_opt(1 Gpc){Θ}/{ρ},

where ρ_opt(1 Gpc) is the binary's optimal signal-to-noise ratio at 1 Gpc.

We additionally simulate cross-correlation measurements of the corresponding stochastic gravitational-wave background, assuming T=1.2 years of integration with Advanced LIGO's Hanford-Livingston baseline.
Our simulated cross-correlation spectra are drawn from
	
    Ĉ(f) ∼𝒩(γ(f)Ω(f), σ(f)),

where the gravitational-wave background's energy density Ω(f) is calculated using Eq. (<ref>) and σ(f) is given by Eq. (<ref>).
Given the binary black hole mass and redshift distributions assumed above and a 1.2 year integration time, the gravitational-wave background has amplitude Ω_0 = 2.2×10^-9 at f=25 Hz and ⟨S/N⟩_opt=4.2.
In our particular noise realization, the binary black hole background is observed with S/N = 3.6, representing a marginal detection.

Figure <ref> illustrates the posterior we obtain on ℛ(z) using our simulated direct detections (left) and direct detections plus stochastic data (right).
Figure <ref> shows the posterior on ℛ_0, α, β, and z_p for this latter case; as before, we have marginalized over the parameters governing the black hole mass distribution.
For reference, Fig. <ref> also includes the one-dimensional marginalized posteriors obtained by direct detections alone (in green).
Full parameter estimation results for each case are given in Table <ref>.

With 500 direct detections we can very precisely measure α=^+_- at 95% credibility, yielding a tight fit to ℛ(z) out to z∼1.
By virtue of not directly observing a turnover of ℛ(z), we can now place a lower limit z_p≥.
Otherwise, we are again limited by Advanced LIGO's finite detection range.
The joint analysis of our direct detections and stochastic data, meanwhile, yields a qualitatively different picture.
Although the S/N of our simulated detection of the gravitational-wave background is somewhat marginal, it provides enough complementary information to rule out large z_p.
While the absolute merger rate remains uncertain at large redshifts, this future data would confidently measure z_p = ^+_-.






§ CONCLUSIONS



We present a powerful new constraint on the binary black hole redshift distribution, with implications for stellar evolution, black hole formation, and binary black hole formation and evolution.
By combining detections of compact binaries in the local Universe with measurements of (or upper limits on) the stochastic gravitational-wave background, we demonstrate that it is possible to explore the binary black hole redshift distribution at redshifts well beyond the present horizon of direct detections.
Using existing observations from the Advanced LIGO/Virgo O1 and O2 observing runs, we have obtained novel joint constraints on the low-redshift slope α and peak z_p of the binary black hole merger rate [see Eq. (<ref>)].
In particular, we can reject merger rates that grow faster than α≳ 7 and peak beyond z_p ≳ 1.5.
These constraints will significantly improve with continued observation.
Given an approximately year-long observation period with design-sensitivity Advanced LIGO, we have demonstrated the possibility of directly measuring z_p.

Although we have taken adopted a decidedly phenomenological model for the merger rate ℛ(z) in this work, this is not the only possible approach.
If, for instance, one were willing to assume that binary black hole formation is tied directly to the (potentially metallicity-dependent) star formation rate, as in Fig. <ref>, one could instead seek to parametrize and measure the metallicity distribution of binary progenitors and the time delay distribution between binary formation and merger.

Looking ahead, future proposed ground-based gravitational-wave detectors like Cosmic Explorer and Voyager may be able to directly measure the rate of binary black hole mergers out to z≳ 10 <cit.>.
However, even a more limited ability to explore the history of binary black hole mergers with present-day instruments will allow us to ask, sooner rather than later, questions of considerable astrophysical importance:
What are the progenitors of compact binary mergers, and when did they form?
What is the mean time delay between binary formation and merger?
How do black hole mergers across cosmic time connect to the evolution of stars and galaxies in the Universe?
The combination of individually resolved sources and the unresolved stochastic gravitational-wave background may soon provide answers.





§ ACKNOWLEDGEMENTS


We would like to thank Nelson Christensen, Andrew Matas, and others within the LIGO Scientific Collaboration and Virgo Collaboration for helpful comments and conversation.
We additionally thank the anonymous referee, whose questions and feedback greatly improved the quality of this work.
TC and WMF thank the Simons Foundation for its generous support of the Flatiron Institute.
TC was also partially supported by the Josephine de Karman Fellowship Trust.
MF was supported by the NSF Graduate Research Fellowship Program under grant DGE-1746045. MF and DEH were supported by  NSF grant PHY-1708081.
They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through an endowment from the Kavli Foundation.
DEH also gratefully acknowledges support from the Marion and Stuart Rice Award.
This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.
LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.



aasjournal.bst


