Does the Black Hole Merger Rate Evolve with Redshift?

The redshift dependence of the binary black hole (BBH) merger rate carries information about the processes by which BBHs evolve and merge, including the environments in which they form, the star formation rate (SFR), and the time-delay distribution. By measuring the luminosity distance⁷ to detected sources, the current generation of ground-based gravitational-wave (GW) detectors will be able to measure the redshift distribution of BBH mergers up to redshifts z ∼ 1 (Abbott et al. 2016b, 2018). The inferred redshift distribution will provide important clues regarding the BBH formation channel. For example, in the classical isolated binary evolution channel, the redshift evolution follows the SFR convolved with a distribution of time delays between formation and merger (Dominik et al. 2015; Abbott et al. 2016b; Belczynski et al. 2016b). Meanwhile, in the dynamical formation channel, the evolution of BBH mergers is tied to the evolution of globular clusters (Chatterjee et al. 2017; Rodriguez et al. 2018). If BBHs are primordial, they are expected to largely follow the dark matter distribution (Mandic et al. 2016; Koushiappas & Loeb 2017). Furthermore, several exotic scenarios, such as gravitational leakage (Dvali et al. 2000; Deffayet & Menou 2007; Pardo et al. 2018) or a significant population of strongly lensed BBH systems (Broadhurst et al. 2018; Smith et al. 2018), would leave an imprint on the inferred redshift distribution.

In this Letter we consider the BBH merger rate density as a function of the component masses, m₁ and m₂, and redshift, z. The mass and redshift distributions must be fit simultaneously, because the detection efficiency of GW detectors depends on the component masses as well as the distance to the source, as discussed in Section 2. We parametrize the mass distribution as a power law with a variable upper mass cutoff, as in Fishbach & Holz (2017). For the redshift distribution, we consider two parametrizations. The first parametrization is motivated by the low-redshift SFR, and assumes that the BBH merger rate follows the comoving volume to zeroth order in redshift. This model can fit astrophysically motivated redshift distributions, including metallicity-weighted SFRs convolved with various time-delay distributions. The second model allows for much more extreme deviations from a uniform in comoving volume merger rate, even at low redshifts, making it less suitable for distinguishing between different formation channels but more sensitive to the exotic scenarios discussed above, including modified gravity. These two models are described in Section 3.

In Section 4.1 we fit a joint mass-redshift distribution to the first six BBH detections announced by the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo Collaboration (LVC): GW150914, LVT151012, GW151226, GW170104, GW170814, and GW170608 (Abbott et al. 2016a, 2016d, 2016e, 2017a, 2017b, 2017c). For simplicity we treat LVT151012, which has an 87% probability of having astrophysical origin, as a full detection (Abbott et al. 2016c). Although it has recently been suggested otherwise (Broadhurst et al. 2018), we find that this set of detections is entirely consistent with a redshift distribution that is uniform in comoving volume. We show how future detections will improve the measurement of the BBH merger rate as a function of component masses and redshift in Section 4.2. We predict that with a few hundred BBH detections by LIGO/Virgo operating at design sensitivity, the measurement of the BBH redshift distribution will be precise enough to distinguish between different formation channels.

The distribution of redshifts among detected BBHs depends on the underlying mass distribution of the black holes (BHs). Assuming that the true BBH merger rate is constant in comoving volume, the redshifts of detected BBHs follow the cumulative probability distributions shown in Figure 1, depending on their component masses. We assume detected BBHs are those that produce a signal-to-noise ratio (S/N) ρ > 8 in a single detector (see Section 4 for more details). Within the relevant mass range for stellar-mass BBHs, the GW signal from more massive BBH mergers is intrinsically stronger ("louder") and can be detected at greater distances. This means that the average redshift among detected heavy BBHs is higher than the average redshift among detected light BBHs. This can be seen in Figure 1, as the cumulative probability curves shift to the right with increasing BBH mass. Equivalently, the mass distribution of detected BHs is different from the true underlying mass distribution, with a preference for more massive BHs over less massive ones. Furthermore, as the sensitivity of the GW detector improves, the average redshift of the detected BBHs will increase. This is seen in the difference between the dashed curves, which assume a noise level appropriate to advanced LIGO's (aLIGO's) second observing run, and the solid curves, which assume the noise level for aLIGO at design sensitivity (respectively, the "Early High Sensitivity" scenario and "Design Sensitivity" noise curves from Abbott et al. 2018).

Zoom In Zoom Out Reset image size

It is clear from Figure 1 that it is impossible to infer the underlying BBH redshift distribution independently of the BBH mass distribution. If the BBH merger rate density (rate per comoving volume) increases with increasing redshift, all of the cumulative probability curves in Figure 1 would shift to the right. However, increasing the relative number of massive BBH mergers in the population also increases the proportion of sources detected at high redshift. In other words, the measured redshift distribution alone cannot distinguish between a merger rate that increases with redshift and a population with a high fraction of massive BBHs (see also Bai et al. 2018). Likewise, the detected mass distribution is sensitive to the underlying redshift distribution; a greater fraction of detected BBHs will be high mass if the merger rate density increases with redshift, because only the high-mass BBHs are detectable at high redshift. In the following, we jointly examine the mass and redshift distribution of merging BBHs.

We consider the differential mass-redshift distribution of BBHs:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}\equiv R\,p({m}_{1},{m}_{2},z),\end{eqnarray} \tag{ 1 }$$

where R is the total number of BBHs across all masses and redshifts, so that

$$\begin{eqnarray}&&\int \displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}{{dm}}_{1}{{dm}}_{2}{dz}=R,\end{eqnarray} \tag{ 2 }$$

and p(m₁,m₂,z) integrates to unity.

Given $\tfrac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}$, we can solve for the usual merger rate density, $\tfrac{{dN}}{{{dV}}_{c}{{dt}}_{m}}$, where t_m is the source-frame time. The merger rate density as a function of redshift is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dV}}_{c}{{dt}}_{m}}(z)=\displaystyle \frac{{dN}}{{dz}}(z){\left[\displaystyle \frac{{{dV}}_{c}}{{dz}}(z)\right]}^{-1}\displaystyle \frac{1+z}{{T}_{\mathrm{obs}}},\end{eqnarray} \tag{ 3 }$$

where:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dz}}=\int \displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}{{dm}}_{1}{{dm}}_{2},\end{eqnarray} \tag{ 4 }$$

T_obs is the total observing time of the GW detector network as measured in the detector frame, and the (1+z) factor converts detector-frame time to source-frame time.

As a first step, we assume that the underlying mass distribution does not vary across cosmic time, so that we can factor the joint mass-redshift distribution as

$$\begin{eqnarray}&&p({m}_{1},{m}_{2},z)=p({m}_{1},{m}_{2})p(z).\end{eqnarray} \tag{ 5 }$$

This assumption may break down over a large range of redshifts, as many formation scenarios predict some dependence of the mass distribution on the merger redshift. However, aLIGO is only sensitive to redshifts z ≲ 1.5 (see Figure 1), where Equation (5) is likely a good approximation, particularly if the distribution of delay times is broad (see, for example, Figure 3 in Chatterjee et al. 2017; Mapelli et al. 2017). For the mass distribution, we use the two-parameter model from Fishbach & Holz (2017), which is an extension to the power-law model employed by the LVC to fit the BBH mass distribution to the first four detections (Abbott et al. 2016c, 2017a) incorporating the possibility of a mass gap above ≳40 M_⊙ due to pair-instability supernovae (Fowler & Hoyle 1964; Heger & Woosley 2002; Belczynski et al. 2016a). We assume that the mass distribution takes the form

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}| \alpha ,{M}_{\max })\propto \displaystyle \frac{{m}_{1}^{-\alpha }}{{m}_{1}-5\ {M}_{\odot }}{ \mathcal H }({M}_{\max }-{m}_{1}),\end{eqnarray} \tag{ 6 }$$

where ${ \mathcal H }$ is the Heaviside step function. We fix the distribution of secondary masses, m₂, to be uniform between the minimum BH mass and m₁, and fix the minimum BH mass, M_min = 5 M_⊙.

3.1. Redshift Model A

For our first redshift model, we choose the following parametrization:

$$\begin{eqnarray}&&p(z| \lambda )\propto \displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{1}{1+z}{(1+z)}^{\lambda },\end{eqnarray} \tag{ 7 }$$

so that λ = 0 reduces to a merger rate density that is uniform in comoving volume and source-frame time. The extra factor of (1+z)⁻¹ converts from detector-frame to source-frame time. Note that for very small z, Equation (7) reduces to a constant in comoving volume and source-frame time merger rate regardless of the value of λ.

If the rate density follows the specific SFR, we would expect

$$\begin{eqnarray}&&p(z)\propto \displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{1}{1+z}\psi (z),\end{eqnarray} \tag{ 8 }$$

where ψ(z) is the specific SFR (Madau & Dickinson 2014)

$$\begin{eqnarray}&&\psi (z)=0.015\displaystyle \frac{{(1+z)}^{2.7}}{1+{[(1+z)/2.9]}^{5.6}}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}.\end{eqnarray} \tag{ 9 }$$

Other models for the SFR, such as Vangioni et al. (2015) or Strigari et al. (2005), agree with the Madau–Dickinson SFR at the low redshifts relevant to aLIGO, z < 1.5. We note that Equation (7) with λ = 2.7 approximates Equation (8) for z ≪ 1, whereas λ = 2.4 provides a very good approximation to Equation (8) for 0.1 ≲ z ≲ 1. Alternatively, because BBH formation is more efficient at low metallicities, we might expect that the rate density follows the low-metallicity SFR, ψ(z)f_Z(z), where f_Z(z) is the fraction of star formation occurring at metallicity ≤Z at redshift z (Belczynski et al. 2010; Mapelli et al. 2013; Spera et al. 2015; Abbott et al. 2016b; Mandel & de Mink 2016). For example, Langer & Norman (2006) give the fit

$$\begin{eqnarray}&&{f}_{Z}(z)=\hat{{\rm{\Gamma }}}\left(0.84,{\displaystyle \frac{Z}{{Z}_{\odot }}}^{2}{10}^{0.3z}\right),\end{eqnarray} \tag{ 10 }$$

where $\hat{{\rm{\Gamma }}}$ is the incomplete gamma function (see also Mandel & de Mink 2016). As the average metallicity decreases with increasing redshift, the low-metallicity SFR rises more steeply with increasing redshift, and peaks at higher redshift. We find that a rate density that follows the low-metallicity (Z ≤ 0.3 Z_⊙) SFR,

$$\begin{eqnarray}&&p(z)\propto \displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{1}{1+z}\psi (z){f}_{Z=0.3{Z}_{\odot }}(z),\end{eqnarray} \tag{ 11 }$$

leads to a redshift distribution, p(z), that is well approximated by Equation (7) with λ = 3.3. It has also been proposed that the progenitors of BBHs are Population III stars formed at zero metallicity, in which case we might expect an even steeper increase of the merger rate with increasing redshift (Belczynski et al. 2004; Kinugawa et al. 2014).

More realistically, the rate density follows the SFR convolved with a time-delay distribution. Different formation channels predict different time-delay distributions. If typical time delays are very long (∼4–11 Gyr), as in the chemically homogeneous formation channel, the rate density will peak at very low redshifts (z ∼ 0.4) well within the aLIGO horizon (Mandel & de Mink 2016). Under the parametrization of Equation (7), this corresponds to λ < 0; in the range z ≤ 1, the best fits to such redshift distributions are given by −6 ≤ λ ≤ −4. In the classical field formation scenario, typical time delays are much shorter (∼10–300 Myr; Dominik et al. 2012, 2013). In this field formation channel, the time delay is expected to follow a distribution

$$\begin{eqnarray}&&\tau \propto {\tau }^{-1}\,{\tau }_{\min }\lt \tau \lt {\tau }_{\max },\end{eqnarray} \tag{ 12 }$$

where typically τ_min ∼ 50 Myr and τ_max is a Hubble time. In the redshift range of interest to aLIGO, this corresponds to a merger density that increases with increasing redshift (λ > 0), but is less steep than the SFR. For example, if the formation rate of BBHs follows the Madau–Dickinson SFR, and the time-delay between formation and merger follows Equation (12), the merger rate at z ≲ 1 can be described by Equation (7) with λ ∼ 1.3. A measurement of λ > 1.3 for the merger rate would indicate that the BBH formation rate density peaks at higher redshift than the SFR (possibly because of metallicity evolution), or that there is a stronger preference for very short time delays.

If the time-delay distribution is in fact restricted to very short time delays (for example, a flat distribution between τ_min = 50 Myr and τ_max = 1 Gyr), the BBH merger rate density will be nearly identical to their formation rate density at z ≲ 1. This may be the case for mergers that take place inside globular clusters. On the other hand, binaries that form dynamically inside clusters but are ejected prior to merger tend to have much longer time delays, on the order of ∼10 Gyr (Rodriguez et al. 2016), corresponding to λ ∼ −10.

3.2. Redshift Model B

For our second redshift model, we assume that the merger rate distribution is uniform in ${V}_{c}^{\gamma /3}={D}_{c}^{\gamma }$ and source-frame time, so that γ = 3 implies that the distribution is uniform in comoving volume. In other words, the redshift distribution takes the form

$$\begin{eqnarray}&&p(z| \gamma )\propto \displaystyle \frac{1}{1+z}\displaystyle \frac{{D}_{c}^{\gamma -1}}{E(z)},\end{eqnarray} \tag{ 13 }$$

where D_c is the comoving distance and ${{dD}}_{c}\propto \tfrac{{dz}}{E(z)}$ (Hogg 1999). This redshift model is more flexible than Model A in the local universe, as it allows for large deviations from a constant in volume merger rate at low redshifts, whereas Model A always reduces to a constant in volume merger rate in the limit z → 0. Model B can constrain scenarios that would cause extreme variations in the slope of the redshift distribution locally, such as a significant population of strongly lensed sources that appear closer than they are (Broadhurst et al. 2018), leading us to infer γ < 3, or GW leakage causing sources to appear farther than they are, leading us to infer γ > 3 (Deffayet & Menou 2007). Previous studies have explored such scenarios through their effects on the S/N distribution (Chen & Holz 2014; Calabrese et al. 2016; García-Bellido et al. 2016) or GW standard siren measurements (Pardo et al. 2018). However, such effects on the redshift distribution (and likewise, the S/N distribution) will most likely be difficult to disentangle from the astrophysical processes that control the redshift and mass distributions. For example, in Figure 2, the dotted pink curve, corresponding to Redshift Model B, γ = 4, and the solid green curve, corresponding to Redshift Model A, λ = 3, are very similar.

Zoom In Zoom Out Reset image size

Figure 2 shows the expected redshift PDF for sources detected by aLIGO at design sensitivity, assuming that the true mass-redshift distribution is described by the models discussed in this section. The solid, dashed, and dashed–dotted blue curves assume the same underlying redshift distribution (corresponding to a constant merger rate density), but different mass distributions. Meanwhile, the solid blue, orange, and green curves show how the detected redshift PDF varies with different underlying redshift distributions parametrized by λ, for a fixed mass distribution. The dotted pink curve assumes the same mass distribution, but takes the underlying redshift distribution to follow Model B with γ = 4. If the merger rate increases with redshift (for example, the solid green compared to the solid blue curve), the detected distribution will skew to high redshifts. However, the effects of changing the mass distribution can be equally, if not more, significant (for example, the difference between the dashed, dashed–dotted and solid blue curves). As we shall see in the following section when we infer the parameters of the mass-redshift model, this leads to a degeneracy between the mass parameters and the redshift evolution parameter.

In this section, we fit our parametrized model for the differential mass-redshift distribution, $\tfrac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}$, to real and simulated LIGO/Virgo detections. Our goal is to extract the four population parameters of the model from GW measurements of the masses and luminosity distances of detected sources. We assume a fixed ΛCDM cosmology determined by the 2015 Planck cosmological parameters (Ade et al. 2016), so that the measured luminosity distance is a direct measurement of the redshift. The shape of the mass-redshift distribution is governed by three parameters, ${\boldsymbol{\theta }}$. For Model A of the redshift evolution, ${\boldsymbol{\theta }}=\{\alpha ,{M}_{\max },\lambda \}$ and for Model B, ${\boldsymbol{\theta }}\,=\{\alpha ,{M}_{\max },\gamma \}$. The fourth parameter, R, corresponds to the total number of detected BBH systems and gives the overall normalization according to Equation (2), allowing us to solve for the physical merger rate of BBHs (Equation (3)).

We model the rate density $\tfrac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}$ as a Poisson point process (Loredo 2004; Farr et al. 2014; Wysocki et al. 2018). The likelihood for the GW data $\{{d}_{i}{\}}_{i=1}^{{N}_{\mathrm{obs}}}$ from N_obs observations, given population parameters $\{{\boldsymbol{\theta }},R\}$ is given by

$$\begin{eqnarray}&&p(\{{d}_{i}\}| {\boldsymbol{\theta }},R)=\left[\displaystyle \prod _{i=1}^{{N}_{\mathrm{obs}}}\int p({d}_{i}| {m}_{1},{m}_{2},z)\displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}\right.\\ &&\quad ({m}_{1},{m}_{2},z| {\boldsymbol{\theta }},R){{dm}}_{1}{{dm}}_{2}{dz}]{e}^{-\beta ({\boldsymbol{\theta }},R)}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&=\left[\displaystyle \prod _{i=1}^{{N}_{\mathrm{obs}}}{\left\langle \displaystyle \frac{1}{\pi ({z}^{i},{m}_{1}^{i},{m}_{2}^{i})}\displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}({m}_{1}^{i},{m}_{2}^{i},{z}^{i}| {\boldsymbol{\theta }},R)\right\rangle }_{\{{z}^{i},{m}_{1}^{i},{m}_{2}^{i}\}}\right.\\ &&\quad \times \,]{e}^{-\beta ({\boldsymbol{\theta }},R)},\end{eqnarray} \tag{ 15 }$$

where $\langle \cdots \rangle$ denotes an average over $\{{z}^{i},{m}_{1}^{i},{m}_{2}^{i}\}$ posterior samples from the ith event, and $\pi ({z}^{i},{m}_{1}^{i},{m}_{2}^{i})$ denotes the interim prior used in the analysis of individual events. The standard priors used in the LIGO analysis of individual events are uniform in component masses and "volumetric" in distance (Veitch et al. 2015):

$$\begin{eqnarray}&&\pi (z,{m}_{1},{m}_{2})\propto {d}_{L}{(z)}^{2}.\end{eqnarray} \tag{ 16 }$$

It is interesting to note that the redshift distribution described by Equation (16) matches Model A with λ = 3 (as opposed to λ = 0, which corresponds to a constant rate density), and so the simplifying assumption that the universe follows a Euclidean geometry implies a redshift distribution that mimics the SFR. Meanwhile, $\beta ({\boldsymbol{\theta }},R)$ is given by

$$\begin{eqnarray}\beta ({\boldsymbol{\theta }},R) & = & \int \displaystyle \frac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}({m}_{1},{m}_{2},z| {\boldsymbol{\theta }},R)\\ & & \times \,{P}_{\det }({m}_{1},{m}_{2},z){{dm}}_{1}{{dm}}_{2}{dz}\end{eqnarray} \tag{ 17 }$$

where P_det(z, m₁,m₂) is the fraction of binary sources at a given redshift and of given component masses that are detectable by the GW detector network. We assume that sources are isotropically distributed on the sky, and the binary inclination is uniformly distributed on the sphere (i.e., uniform in $\cos ({\rm{inclination}})$). We also fix all BH spins to zero in our analysis. Alternatively, we could allow the spins to vary and marginalize over the spin distribution when measuring the mass-redshift distribution (Wysocki et al. 2018). However, incorporating the spin distribution will not affect our analysis significantly, considering that BBHs seem to have small aligned-spin components (Farr et al. 2017, 2018).

Given the component masses and spins, sky position, inclination, and distance (or equivalently, redshift) of the BBH source relative to a GW detector, together with a power spectral density (PSD) that characterizes the noise of the detector, we can calculate the single-detector S/N (Finn & Chernoff 1993; Dominik et al. 2015). We consider a single-detector S/N threshold ρ_th = 8 for detection (corresponding to a network S/N threshold of 12), and assume that for O1 and O2, the noise follows the PSD given by the aLIGO "Early High Sensitivity" scenario (Abbott et al. 2018). For this calculation we ignore the distinction between the true and measured S/N, which, with only six events, does not significantly impact our results. The probability of detection, P_det(z, m₁,m₂), is therefore the fraction of sources that produce a true S/N of ρ = 8 in a single detector.

We are interested in the posterior probability of the population parameters $\{{\boldsymbol{\theta }},R\}$, which is related to the likelihood in Equation (14) by a prior

$$\begin{eqnarray}&&p({\boldsymbol{\theta }},R| \{{d}_{i}\})\propto p(\{{d}_{i}\}| {\boldsymbol{\theta }},R)p({\boldsymbol{\theta }},R).\end{eqnarray} \tag{ 18 }$$

We choose broad, uninformative priors. We take a flat prior for the parameters that make up ${\boldsymbol{\theta }}$ (the power-law slope, maximum component mass, and redshift evolution parameter). Our default prior ranges are α ∈ [−4, 5], M_max ∈ [31 M_⊙, 100M˙_⊙], λ ∈ [−50, 30], and γ ∈ [0, 8]. For the rate parameter, R, we take a flat-in-log prior over the range R ∈ [10,10¹²]. Recall that R is the total number of mergers between redshift z = 0 and the maximum redshift at which the detectors are sensitive—roughly z = 0.6 for O1 and O2 (Abbott et al. 2016b; Kissel et al. 2016). The combined prior is then

$$\begin{eqnarray}&&p({\boldsymbol{\theta }},R)\propto \displaystyle \frac{1}{R}.\end{eqnarray} \tag{ 19 }$$

With this prior choice, the posterior marginalized over R reduces to

$$\begin{eqnarray}&&p({\boldsymbol{\theta }}| \{{d}_{i}\})\\ &&\quad \times \,\propto \displaystyle \prod _{i=1}^{{N}_{\mathrm{obs}}}\displaystyle \frac{\int p({d}_{i}| {m}_{1},{m}_{2},z)p({m}_{1},{m}_{2},z| {\boldsymbol{\theta }}){{dm}}_{1}{{dm}}_{2}{dz}}{\xi ({\boldsymbol{\theta }})}p({\boldsymbol{\theta }}),\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&\xi ({\boldsymbol{\theta }})=\int p({m}_{1},{m}_{2},z| {\boldsymbol{\theta }}){P}_{\det }({m}_{1},{m}_{2},z){{dm}}_{1}{{dm}}_{2}{dz}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&=\,\beta ({\boldsymbol{\theta }},R)/R,\end{eqnarray} \tag{ 22 }$$

and $p({m}_{1},{m}_{2},z| {\boldsymbol{\theta }})$ is related to $\tfrac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}({m}_{1},{m}_{2},z,| {\boldsymbol{\theta }},R)$ by Equation (1). Equation (20) follows because Equation (18) can be written as

$$\begin{eqnarray}&&p({\boldsymbol{\theta }},R| \{{d}_{i}\})\propto p({\boldsymbol{\theta }}| \{{d}_{i}\}){[\xi ({\boldsymbol{\theta }})]}^{{N}_{\mathrm{obs}}}{R}^{{N}_{\mathrm{obs}}-1}{e}^{-R\xi ({\boldsymbol{\theta }})},\end{eqnarray} \tag{ 23 }$$

which when marginalized over R, yields $({N}_{\mathrm{obs}}-1)!\,p({\boldsymbol{\theta }}| \{{d}_{i}\})$. Equation (20) is identical to the form of the posterior derived in Mandel et al. (2016) and used in previous population analyses of GW events (Abbott et al. 2016f; Fishbach & Holz 2017).

4.1. LIGO/Virgo Detections

We fit our mass-redshift model to the first six announced BBH detections. We caution that at the time of writing, the analysis of LIGO's second observing run is still ongoing, and the sample of detections is not guaranteed to be complete. In order to avoid introducing any unmodeled selection biases, proper analysis should wait until the release of the final sample; nevertheless, our analysis illustrates the types of constraints that we expect from six BBHs.

From each of these six events, we approximate the mass and redshift posterior PDFs using the published central values and 90% credible bounds. Specifically, we approximate the detector-frame chirp mass posterior PDFs by Gaussian distributions with a mean and standard deviation that match the published medians and 90% credible widths (Abbott et al. 2016c, 2017a, 2017b, 2017c). The detector-frame masses differ from the source-frame masses by a factor of (1+z) (Krolak & Schutz 1987; Holz & Hughes 2005). Similarly, we use the published medians and 90% credible bounds on the mass ratio, $q=\tfrac{{m}_{2}}{{m}_{1}}$, to find the median and 90% credible bounds on the symmetric mass ratio

$$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{q}{{(1+q)}^{2}}.\end{eqnarray} \tag{ 24 }$$

We then approximate the symmetric mass-ratio posteriors by Gaussian distributions with means and standard deviations that match these medians and 90% credible intervals. We use these approximate posteriors on chirp mass and symmetric mass ratio to approximate the detector-frame component mass posterior distributions for each event. Lastly, we approximate the redshift posterior for each event by a Gaussian distribution matching the published median and 90% credible intervals. We therefore generate posterior samples for the detector-frame masses and redshifts following these approximate distributions. To get posterior samples for the source-frame masses, m₁ and m₂, we divide the posterior samples for the detector-frame masses by (1+z), where the redshifts, z, are drawn from the redshift posterior distribution. This captures the correlations between redshift and source-frame masses in the posterior PDF for an individual event, $p({m}_{1},{m}_{2},z| {d}_{i})$.

Figure 3 shows the resulting posterior PDF on the power-law slope, α, and the redshift evolution parameter, λ or γ, marginalized over M_max and R

$$\begin{eqnarray}&&p(\lambda ,\alpha | \{{d}_{i}\})=\int p(\lambda ,\alpha ,{M}_{\max },R| \{{d}_{i}\}){{dM}}_{\max }{dR}.\end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

Under Model A, we infer $\lambda =-{10}_{-21}^{+15}$, α = 0.7 ± 2.0, and under Model B, we infer $\gamma ={2.7}_{-1.3}^{+1.8}$, $\alpha ={1.1}_{-2.5}^{+2.1}$. All credible intervals are quoted as the median and symmetric (equal-tailed) 90% range. Meanwhile, from the marginal posterior PDF on M_max, we infer that the maximum component BH mass is ${39}_{-6}^{+30}\ {M}_{\odot }$ (Model A) or ${39}_{-6}^{+28}\ {M}_{\odot }$ (Model B); the 95% upper limit of ∼69 M_⊙ is tighter than the 95% upper limit of ∼77 M_⊙ found in Fishbach & Holz (2017) with the additional two detections analyzed here.

There is a positive correlation between the mass power-law slope and the redshift evolution parameter in the two-dimensional posterior, because a mass distribution that favors low masses (large α) is compatible with the data only if the merger rate increases with redshift (large λ or γ) and vice versa; otherwise, more high-redshift and low-mass objects would have been detected. The inferred rate parameter, R, is also positively correlated with these parameters, large α and/or λ (equivalently, γ) imply that there are many more high-redshift low-mass sources that contribute to the total number of mergers, R, but not to the detected number, N_obs (see also Wysocki et al. 2018).

It has recently been suggested that there are statistically too many nearby BBH detections (or equivalently, too many high-S/N detections) compared to the expected constant in comoving volume distribution (Broadhurst et al. 2018). Although our analysis shows a slight preference for a merger rate density that declines with increasing redshift, we find that for both models A and B of the redshift evolution, the current data is consistent with a uniform in comoving volume rate density (λ = 0 or γ = 3) within 1σ (68% credibility). Furthermore, it is possible that the set of published LIGO/Virgo detections is incomplete at this time. If loud events are published first, it is possible that an incomplete set would be biased toward low-redshift events, and our analysis would be artificially biased to small values of λ and γ. A more complete analysis can take place once the results from LIGO's second observing run are finalized.

Finally, we can calculate a posterior PDF on the merger rate density as a function of redshift according to Equation (3). As an illustration of the method, we assume that the total observing time from aLIGO's first and second observing runs is 94 days. This assumption is not based on the true observing time, which is not yet known as analysis on O2 is ongoing. Instead, it is chosen to match the most recently published merger rate estimate from LIGO in Abbott et al. (2017a), which is in the range [12, 213] Gpc⁻³ yr⁻¹ for a BBH population with a uniform in comoving volume merger rate and a mass distribution with power-law slope 1 ≤ α ≤ 2.35. With six published detections, an observing time of 94 days yields a mean merger rate of 100 Gpc⁻³ yr⁻¹ for the "power-law" (α = 2.35) population considered in Abbott et al. (2017a). (Note that this "power-law" mass distribution from Abbott et al. (2017a) fixes the maximum component BH mass to M_max = 95 M_⊙, the minimum component BH mass to M_min = 5 M_⊙, and the maximum total binary mass to M_tot,max = 100 M_⊙.) With this assumption of the observing time, Figure 4 shows the inferred rate density (marginalized over all population parameters) as a function of redshift for Model A (left panel) and Model B (right panel) of the redshift evolution. While previously published merger rate estimates are valid only for a fixed redshift distribution, this method allows us to infer the merger rate simultaneously with the mass and redshift distributions. Once again, we see that our results are consistent with a non-evolving merger rate, which would correspond to a flat horizontal line in Figure 4.

Zoom In Zoom Out Reset image size

From Figure 4, we see that the merger rate is well-constrained at redshifts 0.05 ≲ z ≲ 0.15. Meanwhile, there is too little volume at z ≲ 0.05 to constrain the merger rate well, and the detectors are not sensitive enough at high redshifts. The uncertainties on the merger rate become especially large at high redshifts for Model A, because the parameter for this model, λ, is hard to constrain with low-redshift observations (where the merger rate always approaches a constant in comoving volume rate) but has a significant effect on the high-redshift rate. Meanwhile, varying the parameter, γ, in Model B causes a large variation in the low-redshift rate. This means that γ is easier to measure with low-redshift observations than λ, and so assuming Model B, the merger rate is relatively well-constrained at high redshifts where the sensitivity of the GW detectors approaches zero.

4.2. Future Detections

In this section, we simulate detections from a mock population of BBHs and apply our method to infer the underlying mass and redshift distribution parameters. Our goal is to estimate how many LIGO/Virgo detections will be required to correctly infer a deviation from a uniform in comoving volume merger rate, or alternatively, how many detections will be required to confidently rule out strong deviations from a uniform in comoving volume merger rate. We consider two simulated populations. Both populations follow the same distribution for the BBH masses (Equation (6)), with a minimum component mass of 5 M_⊙, a maximum component mass of 40 M_⊙, and a power-law slope α = 1. We assume that the mass distribution for both populations is independent of the redshift distribution (Equation (5)). The redshift distribution of the first population follows Model A (Equation (7)) with λ = 3, as might be expected if the merger rate followed the low-metallicity SFR convolved with a time-delay distribution that favored short time delays. Meanwhile, the second population has a uniform in comoving volume merger rate, corresponding to λ = 0 in Model A or γ = 3 in Model B.

For each mock population we generate BBHs with component masses and redshifts following the assigned underlying distribution. Only a subset of BBHs are detected, and their measured masses and redshifts take the form of (marginalized) posterior PDFs. In order to generate realistic mass and redshift measurements, we construct a synthetic detection model. The synthetic detection model enables us to self-consistently and realistically capture the correlations between the measured S/N, which determines the detectability of an event and its measured redshift. We verify that the model closely approximates the detectability of a BBH with given masses and redshift by aLIGO at design sensitivity.

Synthetic BBHs are generated as follows. Each BBH system is characterized by four parameters: the source-frame chirp mass ${ \mathcal M }$, the symmetric mass ratio η, the luminosity distance d_L, and an angular factor Θ. Each system also has an associated true S/N, ρ, which depends on these four parameters. Note that ${ \mathcal M }$ and η allow us to directly infer m₁ and m₂, and d_L allows us to infer z. Here, Θ plays the combined role of the sky location, inclination, and polarization on the measured GW amplitude. We tune the width of the Θ distribution to control the uncertainty of the measured signal strength, which in turn controls the uncertainty on the measured luminosity distance. This allows us to capture the correlations between the measured signal strength and the measured redshift, which is necessary in order to model the selection effects consistently between the detection model and the calculation of P_det (Equation (17)).

We set the "typical" S/N, ρ₀, of a BBH system with parameters ${ \mathcal M }$ and d_L to

$$\begin{eqnarray}&&{\rho }_{0}\equiv 8{\left(\displaystyle \frac{{ \mathcal M }(1+z)}{{{ \mathcal M }}_{8}}\right)}^{5/6}\displaystyle \frac{{d}_{L,8}}{{d}_{L}},\end{eqnarray} \tag{ 26 }$$

where we fix ${{ \mathcal M }}_{8}=10\ {M}_{\odot }$ and d_L,8 = 1 Gpc. This scaling approximates the amplitude of an inspiral GW signal to first order, and we chose ${{ \mathcal M }}_{8}$ and d_L,8 to roughly match the typical distances of detected sources by aLIGO at design sensitivity (Chen et al. 2017).⁸ The true S/N, ρ, in our model depends on the angular factor Θ, and is given by

$$\begin{eqnarray}&&\rho ={\rho }_{0}{\rm{\Theta }},\end{eqnarray} \tag{ 27 }$$

similar to relationship between the true S/N and the "optimal S/N" via the projection factor Θ (Finn & Chernoff 1993) or w (Dominik et al. 2015; Chen et al. 2017), although Θ in our case is simply a random variable with a log-normal distribution. As the uncertainty on Θ controls the uncertainty on d_L, we pick the variability of Θ in order to get realistic measurement uncertainties on d_L. We find that realistic measurement uncertainties on d_L are achieved when Θ has a typical width of 15%, and so we pick

$$\begin{eqnarray}&&\mathrm{log}{\rm{\Theta }}\sim N\left(0,\displaystyle \frac{0.3}{1+\tfrac{{\rho }_{0}}{8}}\right).\end{eqnarray} \tag{ 28 }$$

From the true parameters ${ \mathcal M }$, η, d_L, and Θ, we assume that the measurement process measures three parameters: the observed S/N, ρ_obs, the observed chirp mass, ${{ \mathcal M }}_{\mathrm{obs}}$, and the observed symmetric mass ratio, η_obs. These are given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{obs}}=\rho +N(0,1),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\mathrm{log}{{ \mathcal M }}_{\mathrm{obs}}=\mathrm{log}({ \mathcal M }(1+z))+N\left(0,8\displaystyle \frac{{\sigma }_{{ \mathcal M }}}{{\rho }_{\mathrm{obs}}}\right),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\eta }_{\mathrm{obs}}=\eta +N\left(0,8\displaystyle \frac{{\sigma }_{\eta }}{{\rho }_{\mathrm{obs}}}\right),\end{eqnarray} \tag{ 31 }$$

where we assume that the observed S/N is normally distributed about the true S/N, ρ, with a standard deviation of 1 (due to different realizations of Gaussian noise), and the uncertainties on the mass parameters scale inversely with the observed S/N (Veitch et al. 2015). We fix a threshold of ρ_obs ≥ 8 for detection. To match the expected measurement uncertainties, we fix ${\sigma }_{{ \mathcal M }}=0.04$ and σ_η = 0.03, so that the relative 90% credible interval uncertainty for the recovered detector-frame primary (secondary) mass is typically 40% (50%; Ghosh et al. 2016; Vitale et al. 2017). Meanwhile, the luminosity distance is measured from ρ_obs via Equations (26)–(28). The typical relative 90% confidence interval uncertainty for the recovered luminosity distance is ∼50%, which is also a realistic expectation (Vitale et al. 2017). As discussed earlier in this section, this process of recovering the measured luminosity distance from the measured signal strength is necessary in order to incorporate selection effects consistently, and ensure that we generate single-event posteriors, p(m₁, m₂, z) that are compatible with the assumed detection probability, P_det(m₁, m₂, z).

Under this synthetic (yet realistic) detection model, we generate 500 detected BBH systems for each of the two populations. The projected constraints on the power-law slope and redshift evolution parameter are shown in Figure 5. Note that the maximum mass parameter will be already tightly measured with a few tens of detections (Fishbach & Holz 2017). We find that after 100 detections by LIGO/Virgo (which may happen as early as the next observing run, starting in late 2018) it may be possible to detect deviations from a uniform in comoving volume merger rate if the true redshift distribution evolves as steeply as the low-metallicity SFR (λ ∼ 3). (This is expected from formation channels where the typical time delay between formation and merger is short.) Additionally, we expect to distinguish between a merger rate density that increases with increasing redshift and a merger rate density that decreases with redshift (as expected from formation channels with very long time delays). If the true deviation from a uniform merger rate density is small (0 < λ < 1, where λ = 0 implies a uniform merger rate density), it may take ∼500 detections to confidently exclude λ = 0. This will require a few years of aLIGO operating at design sensitivity (starting in 2020+; Abbott et al. 2018). Meanwhile, extreme deviations from a constant merger rate density ($\gamma \ne 3$ in Redshift Model B) can be ruled out in 100 detections, as γ will be constrained to a 90% credible interval width of ≲1.

Zoom In Zoom Out Reset image size

In summary, we expect that with N detections by LIGO/Virgo operating at design sensitivity, we will be able to constrain λ from Model A to a 90% credible interval of width $\sim 31/\sqrt{N}$, and γ from Model B to a 90% credible interval of width $\sim 8.5/\sqrt{N}$, although the exact rate of convergence depends on the true values of these parameters. If the local BBH merger rate is 100 Gpc⁻³ yr⁻¹, the mass distribution follows a power law with α = 1 and M_max = 40 M_⊙, and the redshift distribution is constant in comoving volume in source-frame time, we expect ∼300 detections in one year of LIGO/Virgo operating at design sensitivity (assuming that each detector has a duty cycle of 80%, and that a confident detection requires at least two detectors in observing mode). Fixing this mass distribution, but assuming instead that the redshift distribution follows the Madau–Dickinson SFR, the expected number of detections increases to ∼680. On the other hand, if the overall merger rate is the same, but the BBH mass distribution follows a steeper power law with α = 2.35 (fixing M_max = 40 M_⊙), we expect ∼140 detections in a year of LIGO/Virgo operating at design sensitivity if the merger rate density is independent of redshift, and ∼280 detections if the merger rate density follows the SFR. (For O3 sensitivity, the expected number of detections is smaller by a factor of 4–5.) With 100–700 detections per year, we expect to detect deviations from a constant merger rate density, or severely constrain such deviations, within the first 1–3 years of LIGO/Virgo operating at design sensitivity (starting ∼2020).

We have explored the ability of the LIGO/Virgo network to measure the redshift evolution of the BBH population. We have applied a simple four-parameter model to constrain the BBH merger rate density, $\tfrac{{dN}}{{{dm}}_{1}{{dm}}_{2}{dz}}$, as a function of component masses and redshift. Our model allows us to simultaneously constrain the slope and maximum mass of the distribution of primary BH masses, the slope of the redshift distribution, and the merger rate. We note that our method can also be applied to the binary neutron star (BNS) population, although such sources will only be detectable up to redshifts z < 0.1 with the current generation of GW detectors. However, the mass distribution of such sources may already be well constrained from the galactic population (Özel et al. 2012). If we adopt this mass distribution as a prior, the analysis on BNS would simplify to a one-parameter redshift evolution model.

Recall that a measurement of the merger redshift distribution constrains a combination of the formation rate as a function of redshift and the time-delay distribution. If we can constrain the redshift evolution parameter to λ ≳ 2.4, we may infer that the BBH formation rate density peaks at higher redshift than the SFR, regardless of the time-delay distribution. If we assume that the time-delay distribution follows Equation (12) with τ_min = 50 Myr and τ_max = 14 Gyr, a measurement of λ ≳ 1.3 implies that the BBH formation rate density peaks at higher redshift than the SFR. Alternatively, if we assume that the formation rate density follows the low-metallicity SFR in Equation (11), measuring λ ≳ 1.9 (λ ≲ 1.9) would allow us to infer that the time-delay distribution is more skewed toward short (long) time delays than Equation (12).

Because our focus is on extracting the redshift evolution of the merger rate, we have simplified our treatment of the mass distribution. For example, our parametrization assumes that the BBH mass distribution does not evolve with redshift, and does not allow the distribution of mass-ratios or the minimum BH mass to vary. We have verified that adding two free parameters, β and M_min, to describe the mass-ratio distribution and the minimum mass according to

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}| \alpha ,\beta ,{M}_{\max },{M}_{\min })\propto \displaystyle \frac{{m}_{1}^{\alpha }{m}_{2}^{\beta }}{{m}_{1}-{M}_{\min }}{ \mathcal H }({M}_{\max }-{m}_{1}),\end{eqnarray} \tag{ 32 }$$

does not significantly affect our results for the six BBH detections. Taking uniform priors on these additional parameters, −4 < β < 4, 3 < M_min < 9 M_⊙,⁹ causes the posteriors on the remaining four parameters to widen only slightly. Furthermore, the constraints on the additional two parameters, β and M_min, are not informative with only six detections, and are consistent with our default values, β = 0, M_min = 5 M_⊙, although we find a slight preference for equal mass ratios (β > 0 at ∼70% credibility; consistent with the results of Roulet & Zaldarriaga 2018) and a larger minimum mass (M_min > 3.9 M_⊙ at 95% credibility). These results hold for both Models A and B of the redshift evolution. It will likely take ${ \mathcal O }(100)$ detections to measure M_min sufficiently well and resolve the putative gap between the neutron star and BH mass spectrum (Littenberg et al. 2015; Kovetz et al. 2017; Mandel et al. 2017).

With sufficient detections, our four-parameter model will likely break down, and we should include more degrees of freedom in the mass-redshift (and possibly also spin) distribution to avoid introducing systematic biases in the inferred parameters (Talbot & Thrane 2018; Wysocki et al. 2018). For example, if the mass distribution varies with redshift, possibly favoring larger masses at high redshifts due to the lower average metallicity, our simple model will misinterpret this as an evolution in the merger rate. Therefore, a more complicated model should allow for correlations between the mass and redshift distribution, either through the addition of one to two parameters (e.g., a copula model), or a many-parameter model that fits the mass distribution separately in different redshift bins. It may also be possible to introduce a multi-component mixture model to determine whether there are multiple populations of BBHs following different mass-redshift distributions. Furthermore, a more sophisticated model would include at least two additional parameters to fit the peak, z_peak, and the high-redshift slope of the merger rate density. For example, we can consider the following parametrization of the merger rate density inspired by the Madau–Dickinson SFR:

$$\begin{eqnarray}&&p(z| a,b,{z}_{\mathrm{peak}})\propto \displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{{(1+z)}^{a}}{1+\left(\tfrac{a}{b-a}\right){\left[\tfrac{1+z}{1+{z}_{\mathrm{peak}}}\right]}^{b}}.\end{eqnarray} \tag{ 33 }$$

This parametrization provides an excellent fit to all of the astrophysical redshift distributions discussed in Section 3 up to redshifts z ∼ 4, whereas our one-parameter model of Equation (7) starts to break down at z ∼ 1. However, it is unlikely that we will have tight constraints on z_peak and b with second-generation GW detectors, because z_peak will likely lie beyond the sensitivity of these detectors, unless the time delays between formation and merger are typically extremely long (greater than a few Gyr). For example, if the BBH formation rate follows the Madau–Dickinson SFR and the time-delay distribution in Equation (12), the peak of the merger rate density would be at z_peak = 1.4, where we expect to have very few detections (see Figures 1 and 2). If the BBH formation rate peaks later than the Madau–Dickinson SFR, following the low-metallicity SFR in Equation (11) for example, the peak of the merger rate density would be even farther out of reach, at z_peak = 2.1, assuming the same time-delay distribution. However, by the time we have ${ \mathcal O }(1000)$ detections, it may be possible to get some measurement of z_peak, which would allow us to infer the peak of the formation rate density for an assumed time-delay distribution (this peak would be at z = 1.8 if BBH formation followed the Madau–Dickinson SFR, and z = 2.7 if the formation followed the low-metallicity SFR of Equation (11)). We look forward to a time where our single redshift-evolution parameter is sufficiently well measured that we must include these additional parameters in our model; we anticipate that this will take over 500 detections. We note that with third-generation GW detectors, it will be possible to accurately infer the entire formation rate history of BBHs together with the time-delay distribution from the observed redshift evolution of the merger rate (Vitale & Farr 2018).

In addition to the limitations of our parametrized model, another much less significant source of systematic uncertainty in our analysis comes from GW measurements of the luminosity distance (and therefore, the redshift) to a source. Extracting the luminosity distance from a GW signal depends on measuring its amplitude, which is affected by detector calibration uncertainties. The calibration uncertainty is only a few percent (Karki et al. 2016), which is negligible compared to the expected uncertainty on the redshift evolution parameter. Another subdominant source of uncertainty comes from the effect of weak lensing on the GW amplitude, which contributes at the sub-percent level and is therefore negligible for our analysis (Holz & Wald 1998; Holz & Linder 2005).

By fitting a four-parameter mass-redshift distribution to the first six announced BBH mergers, we have placed the first constraints on the redshift evolution of the BBH merger rate. We show that because of strong correlations between the masses and redshifts of detected BBHs, the mass and redshift distribution must be fit simultaneously. We consider two parametrizations of the redshift evolution: Model A fixes the slope of the redshift distribution to match the differential comoving volume locally (as z → 0), as is expected from most astrophysical formation channels, while Model B allows for large deviations, even at low redshift. Our constraints from six events are too weak to distinguish between any astrophysical formation scenarios: for example, we measure −31 < λ < 5 at 90% credibility for Model A, whereas typical formation channels predict between −4 ≲ λ ≲ 3). However, we can already constrain extreme deviations from a uniform in comoving volume merger rate, finding $\gamma ={2.7}_{-1.3}^{+1.8}$ for Model B (γ = 3 corresponds to a uniform merger rate). Furthermore, while previous analyses have calculated the BBH merger rate under the assumption of a uniform in comoving volume redshift distribution, we demonstrate how to infer the merger rate density as a function of redshift.

We project that with 100–500 detections by LIGO/Virgo, the inferred redshift evolution of the BBH merger rate will allow us to distinguish between proposed formation channels (for example, those that favor long versus short time delays between progenitor formation and BBH merger). Meanwhile, the overall merger rate and mass distribution will also provide important clues regarding the formation channel (Dominik et al. 2015; Stevenson et al. 2015; Zevin et al. 2017; Barrett et al. 2018). In 5–10 years, the constraints on the redshift evolution parameter alone will allow us to infer the peak of the formation rate of BBH progenitors and/or the typical time delay between formation and merger.

We thank Christopher Berry, Eve Chase, Reed Essick, and Peter Shawhan for their helpful comments on the manuscript. M.F. was supported by the NSF Graduate Research Fellowship Program under grant DGE-1746045. M.F. and D.E.H. were supported by NSF grant PHY-1708081. They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through NSF grant PHY-1125897 and an endowment from the Kavli Foundation. D.E.H. also gratefully acknowledges support from the Marion and Stuart Rice Award. W.M.F. is supported in part by the STFC.