Probing Multiple Populations of Compact Binaries with Third-generation Gravitational-wave Detectors

Advanced gravitational-wave (GW) detectors such as LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2015), and Kagra (Aso et al. 2013) have dramatically increased our ability to study stellar mass black holes (BHs) and the environments in which they form. The latest catalog released by the LIGO-Virgo-Kagra (LVK) collaboration includes 39 new GW detections, most of which are binary BHs (BBHs), bringing the total number of stellar mass BHs detected with GWs to over 100 (Abbott et al. 2019a, 2020a). These observations have already allowed for interesting astrophysical measurements (Abbott et al. 2019a, 2019b, 2020a, 2020b), such as hints of multiple formation channels (Farr et al. 2017; Zevin et al. 2017; Antonini & Gieles 2020; Belczynski et al. 2020; Zevin et al. 2021; Callister et al. 2020a; Wong et al. 2021b), hierarchical mergers (Fishbach et al. 2017; Gerosa & Berti 2017; Doctor et al. 2019; Kimball et al. 2019, 2020b; Gerosa et al. 2020; Tiwari & Fairhurst 2020), and constraints on stellar physics (Farmer et al. 2020; Fragione & Loeb 2021; Bavera et al. 2021), primordial BHs (Bird et al. 2016; Sasaki et al. 2016; Ali-Haïmoud et al. 2017; Boehm et al. 2020; Hall et al. 2020; Hütsi et al. 2021; Wong et al. 2021a), and ultralight bosons (Arvanitaki et al. 2017; Brito et al. 2017; Ng et al. 2021a, 2021b) Therefore, the growing set of BBHs enables the study of the properties of both individual sources and the underlying populations.

Among the key questions that can be addressed by GW astrophysics are how many such populations exist and what their characteristics are. Multiple approaches have been proposed to address these questions, which ultimately rely on looking for features that would be expected in the BBH generated by each channel, such as their mass and spin distribution (Stevenson et al. 2015; Vitale et al. 2017; Fishbach & Holz 2017; Farr et al. 2017, 2018; Talbot & Thrane 2017, 2018; Bouffanais et al. 2019; Doctor et al. 2019; Fishbach et al. 2020; Fishbach & Holz 2020; Miller et al. 2020; Safarzadeh et al. 2020; Kimball et al. 2020a), eccentricity distribution (Lower et al. 2018), or redshift distribution (Fishbach et al. 2018; Callister et al. 2020b). However, due to the limited sensitivity of current GW detectors, the observed BBHs are relatively "local," with redshift z ≲ 1 (Abbott et al. 2019a, 2019b, 2020a, 2020b; Roulet et al. 2020; Venumadhav et al. 2020). Even as LIGO, Virgo, and Kagra improve their sensitivities with the implementation of frequency-dependent quantum squeezing (McCuller et al. 2020) and better low-frequency isolation (Yu et al. 2018), the detector horizon will reach z ∼ 3 only for the heaviest systems (Abbott et al. 2016; Hall & Evans 2019).

Second-generation detectors will therefore be unable to access mergers at redshifts larger than a few. This likely precludes the possibility of detecting binaries whose component BHs originated directly from the Population III (Pop III) stars, whose formation peak could be as high as z ∼ 10. Primordial BHs (Kinugawa et al. 2014, 2016; Hartwig et al. 2016; Belczynski et al. 2017; Raidal et al. 2017, 2019) might also be out of reach for existing facilities. ⁶ This will not be the case in the era of third-generation (3G) GW detectors such as the Cosmic Explorer (CE; Abbott et al. 2017; Reitze et al. 2019) and Einstein Telescope (ET; Punturo et al. 2010; Van Den Broeck 2012; Maggiore et al. 2020). In fact, 3G detectors will have horizons up to z ≳ 30 and enable access to most of the BBHs throughout the cosmic history (Hall & Evans 2019). Therefore, both detections and nondetections of high-redshift BBHs with 3G detectors can provide significant constraints on the properties of Pop III remnants. This is particularly important, as the remnants of Pop III stars might be the light seeds that lead to the formation of supermassive BHs early in cosmic history (Greene et al. 2020). Whereas the space-based GW detector LISA (Amaro-Seoane et al. 2017) and electromagnetic missions such as the X-ray observatory Lynx (Gaskin et al. 2019) might reveal the presence of heavy (i.e., M ≳ 100 M_⊙) BHs at redshifts of 10, ground-based 3G detectors might very well be the only way to search for smaller building blocks: stellar mass BHs.

Using the current local BBH merger rate estimate, ∼25 Gpc⁻³ yr⁻¹ (Abbott et al. 2020b), the total merger rate of the BBHs in the universe is inferred to be ∼10,000 month^–1, if the merger rate density follows the same evolution as the star formation rate (SFR; Regimbau et al. 2017). Nearly all of these BBHs would be detectable by 3G detectors, most of which will also have a very high signal-to-noise ratio (S/N), leading to very precise distance (and hence redshift, assuming a known cosmology) measurements (Vitale 2016; Vitale & Evans 2017; Vitale & Whittle 2018). The large number of loud observations thus allows for inferring the morphology of merger rate densities in both parametric and nonparametric ways (Van Den Broeck 2014; Safarzadeh et al. 2019; Vitale et al. 2019; Romero-Shaw et al. 2020). In Vitale et al. (2019), we showed how, by combining redshift measurements of tens of thousands of GW observations in 3G detectors and assuming that all BBHs are formed in galactic fields, one can infer the SFR history and time-delay distribution. In light of the properties of the BBHs detected by LIGO and Virgo in the last few years, it has become harder to assume that all BHs are formed in galactic fields, since many of the BHs being detected are consistent with having formed dynamically, in globular or nuclear clusters, or in the disks of active galactic nuclei (AGN; Bartos et al. 2017; Yi & Cheng 2019; Yang et al. 2019, 2020; Gröbner et al. 2020; Samsing et al. 2020; Tagawa et al. 2020a, 2020b, 2021). In this paper, we greatly extend our previous work (Vitale et al. 2019) and explore the possibility of simultaneously identifying and constraining the properties of multiple formation channels. In particular, we simulate universes where BHs can be formed in galactic fields and globular clusters and from Pop III stars. We show how well one can constrain the properties of each channel and their branching ratios. In particular, we focus on the evidence for a high-redshift (z ≳ 5) population, which would be the smoking gun of formation outside of the traditional evolutionary pathways.

In this section, we summarize the main astrophysical properties of the BBH formation channels that we consider in the paper. Astrophysical BBHs are believed to form in various environments, such as binary stellar evolution in galactic fields (Dominik et al. 2012, 2013, 2015; de Mink & Belczynski 2015; Belczynski et al. 2016; O'Shaughnessy et al. 2017; Stevenson et al. 2017; Broekgaarden et al. 2019; Mapelli et al. 2019; Breivik et al. 2020; Bavera et al. 2020), dynamical formation through multibody interactions in star clusters (from low-mass stars to large nuclear star cluster; Portegies Zwart & McMillan 2000; Rodriguez et al. 2015, 2016; Rodriguez & Loeb 2018; Di Carlo et al. 2019; Antonini & Gieles 2020; Santoliquido et al. 2021; Kremer et al. 2020) or AGN disks (Bartos et al. 2017; Yi & Cheng 2019; Yang et al. 2019, 2020; Gröbner et al. 2020; Samsing et al. 2020; Tagawa et al. 2020a, 2020b, 2021), Population III (Pop III) stars (Kinugawa et al. 2014, 2016; Hartwig et al. 2016; Belczynski et al. 2017), or primordial BHs (Carr & Hawking 1974; Bird et al. 2016; Sasaki et al. 2016; Ali-Haïmoud et al. 2017; Clesse & García-Bellido 2017; Raidal et al. 2017, 2019; Boehm et al. 2020; Hall et al. 2020; Wong et al. 2021a).

For simplicity, in this paper, we use galactic fields and globular clusters as the only main populations at low redshift and Pop III stars as the only high-redshift channel. The analysis can be easily extended to even more populations at the price of an increased computational cost.

Hereafter, we label the BBHs in the three formation channels as field binaries, cluster binaries, and Pop III binaries. Since an astrophysical BH is a remnant of stellar collapse, the merger rate history of each channel is correlated with the SFR and the time delay from binary formation to merger. Cluster and field binaries consist of BH remnants left over from Population I/II stars, whose corresponding SFR peaks at late times: z ∼ 3 (Madau & Dickinson 2014; Vangioni et al. 2015). Accounting for the typical time delay between binary formation and merger (∼10 Myr to ∼10 Gyr), their merger rates are expected to peak at around z ∼ 2 (Dominik et al. 2013; Belczynski et al. 2016; Rodriguez & Loeb 2018; Mapelli et al. 2019).

Pop III stars are instead formed at early times, z ≳ 10, from primordial gas clouds at extremely low metallicity (de Souza et al. 2011; Vangioni et al. 2015). The "metal-free" environment reduces the stellar wind mass loss during the binary evolution (Baraffe et al. 2001) so that Pop III stars might be more massive than later stellar populations. Eventually, heavy BH remnants are left behind that merge in a short timescale, resulting in a merger rate density that could peak at around z ∼ 10 (Kinugawa et al. 2014, 2016; Belczynski et al. 2016; Hartwig et al. 2016).

The fact that different formation channels result in different merger rate distributions as a function of redshift, especially for Pop III remnants, can be exploited to infer the properties of BBH populations solely based on the redshift distribution of the detected sources. More elaborate tests can be envisaged that also rely on other distinguishing features, e.g., masses, spins, or eccentricity of the sources. In this study, we will show that tests based on the redshift distribution alone can already provide significant constraints, while also having the benefit of being model-independent, at least in some of the implementations we demonstrate below. This seems particularly desirable, since the true distribution of intrinsic parameters such as masses and spins is highly uncertain, especially for Pop III remnants.

In Figure 1, we show the "true" merger rate densities of three formation channels we use in this study and focus on key features of their shapes (we will discuss the branching ratios, i.e., the relative scale, later). The rate densities of field (orange), cluster (blue), and Pop III (green) binaries are phenomenological fits (details in Appendix C, Equations (C1)–(C3)) of the population synthesis simulation from Belczynski et al. (2016), Rodriguez & Loeb (2018), and Belczynski et al. (2017), respectively. We notice that the high-redshift tail of the cluster merger rate density is much steeper than that of the field binaries. This is largely due to the choice of model; the cluster merger rates are based on the model of globular cluster formation from El-Badry et al. (2019), which goes to zero at z ∼ 10. That, combined with the delay between cluster formation and BBH mergers (since BHs in clusters can only merge after the cluster has formed and the BHs have sunk to the center due to dynamical friction, a process that can take ∼100 Myr; e.g., Morscher et al. 2015), causes a steeper slope in the merger rate at high z. The field and cluster merger rate densities peak at similar values, z ∼ 2.2 and 2.6, respectively, whereas the Pop III merger rate density peaks much later, at z ∼ 11.6.

Zoom In Zoom Out Reset image size

The vertical lines in Figure 1 report the horizon of future ground-based detector networks. Advanced detector networks can observe BBHs up to the low-redshift peak of the merger rate densities of the two dominating channels, fields, and clusters. However, as the plot shows, for z ≲ 2, the merger rate densities of both field and cluster binaries are quite similar. Hence, advanced detectors are unlikely to be able to disentangle the two channels using only redshift information (as mentioned above, one can use other features, at the price of making the analysis more model-dependent). The situation improves with a single Voyager detector ("Voyager-1" line), which can access most of the field and cluster binaries up to z ≲ 8 and therefore exploit the expected difference in their merger rate after the peak to characterize the two channels. A network of three Voyager-like detectors ("Voyager-3" line) can extend the horizon to a redshift where the contribution to the total merger rate of the field and the Pop III channel might become comparable. However, it is only with 3G detectors that one can access the peak of the merger rate from Pop III. In fact, the horizon of CE and ET to heavy BBHs is outside of the range of Figure 1 (as indicated by the cyan arrow in the bottom right corner), at z ∼ 100. As we will show in the following sections, the fact that the horizon of 3G detectors extends well beyond the expected peak of Pop III mergers allows for both modeled and unmodeled tests of the existence of such a subpopulation.

In this work, we follow two approaches to measure the comoving-frame merger rate density dR/dz: (i) a unmodeled approach that utilizes Gaussian process regression (GPR) to infer dR/dz as a piecewise function over several redshift bins (Mandel et al. 2017) and (ii) a modeled approach in which we use phenomenological models for the various subpopulations (Farr et al. 2015; Vitale et al. 2019). In both cases, we use a hierarchical Bayesian inference framework (Farr et al. 2015; Mandel et al. 2019; Thrane & Talbot 2019; Wysocki et al. 2019; Vitale 2020) to measure the parameters of the population(s). More details are provided in Appendix A. Details about the implementation of the GPR analysis can be found in Appendix B, whereas Appendix C reports the functional forms of the modeled populations. The priors used in the analysis are documented in Appendix E.

To study how well the models can identify the Pop III subpopulation, we perform a mock-data challenge by simulating 18 different universes, which contain 2 months worth of BBH data with a majority of cluster or field binaries. The detailed setup of the simulations can be found in Appendix D.

In the following, we will focus on the measurement of the volumetric merger rate density, $\dot{n}(z)\equiv {dR}/{{dV}}_{c}$, rather than R itself. (We will use an index to indicate the volumetric merger rate in a specific channel, e.g., ${\dot{n}}_{\mathrm{III}}(z)$ for the volumetric merger rate in the Pop III channel. Here "F" will indicate the field channel and "G" the globular cluster channel.) We will also report the branching ratios between the channels. Here f_III will indicate the fraction of Pop III mergers over the total, whereas ${\tilde{f}}_{G}$ will indicate the fraction of cluster binaries over the sum of field and cluster binaries.

Our modeled approach naturally provides more information about the characteristic parameters of each channel and their correlations. Those are discussed in Appendix F.

3.1. Unmodeled Analysis

Figure 2 shows our inference on the volumetric merger rate $\dot{n}$ obtained with the unmodeled GPR approach. The gray bands report the 68% and 95% credible intervals for the universes with (right panel) and without (left panel) Pop III binaries. The colored lines show the true volumetric merger rate of the individual populations. For both of the panels, the true branching ratio between field and cluster binaries is ${\tilde{f}}_{{\rm{G}}}=0.5$. We find that the true $\dot{n}$ values (black solid lines) lie within the 95% credible intervals in both cases. While the relative uncertainty on $\dot{n}$ is at a percent level at z ≲ 6, it increases to ${ \mathcal O }(100 \% )$ at z ≳ 8. This is because (i) the S/N of each source decreases with the distance and (ii) the number of sources in each redshift bin is decreasing as the differential comoving volume shrinks at earlier times in the history of the universe.

Zoom In Zoom Out Reset image size

Perhaps the most attractive feature of the unmodeled analysis is that we can find some evidence for the presence of a high-redshift subpopulation, even without strong modeling. The simplest way of doing this is to look for local peak(s) in $\dot{n}$. At the very minimum, we would expect to find evidence for the "main" peak at z ∼ 2 arising from the merger in fields and clusters, while high-redshift peaks would be indicative of a different subpopulation.

We implement a peak finder algorithm simply by asking the first derivative of $\dot{n}$ to be zero and the second to be negative: $d\dot{n}/{dz}=0$ and ${d}^{2}\dot{n}/{{dz}}^{2}\lt 0$. Some care is required to avoid false positives due to natural oscillations in the results of the GRP that are not due to astrophysical maxima (or minima) but only to the underlying Gaussian process. These are particularly visible at high redshifts in Figure 2(a). To mitigate the effect of these fluctuation, we require the height of any high-redshift peaks to be at least 1/10 of the height of the low-redshift peak, as well as an intrapeak separation larger than Δz = 1. These requirements are based on the expected excess of Pop III, as discussed in Appendix D, and arguably represent the only modeling involved in the GRP approach that we describe.

With these two restrictions, we count the number of peaks N_peak for each $\dot{n}$ sample of the GRP in every simulated universe. The results are shown in Figure 3 as a function of the true branching ratio cluster/field. For the universes without Pop III binaries, we recover a single peak, as expected, with >99% probability (purple histograms). This implies the nonexistence of a secondary peak whose relative height is ≥10% of the primary peak. On the other hand, the true value for the number of peaks, N_peak = 2, is found at ≳90% credibility for the universes with Pop III binaries (yellow histograms).

Zoom In Zoom Out Reset image size

We observe that in the universes with Pop III binaries, the posterior on the number of peaks has a secondary mode at N_peak = 3 for ${\tilde{f}}_{{\rm{G}}}\lesssim 0.5$. We explain this as follows: for ${\tilde{f}}_{{\rm{G}}}\lesssim 0.5$, the field binaries are the dominating channel and, as is clear in Figure 2, produce a flatter high-redshift tail than the cluster channel. Hence, it is easier to produce multiple peaks due to fluctuation and induce a leakage to N_peak ≥ 3. On the other hand, when ${\tilde{f}}_{{\rm{G}}}\gtrsim 0.5$, the cluster binaries are dominating and do not contribute to the merger rate at z ≳ 8, where Pop III becomes the only source of BBHs, which makes the high-redshift peak narrower and hence easier to reveal. But the Poisson fluctuation of the high-redshift bins in dR/dz may still lead to an underestimation of the relative height below our 10% threshold, inducing a small contamination at N_peak = 1. We note that the above trend is subject to the model uncertainty of our chosen simulation data in the high-redshift region.

This method also allows for measuring the location(s) of the peak(s), which would be useful to understand the population properties. For instance, the shift of the primary peak relative to the SFR could inform the typical time delay to merger and a hint of metallicity evolution (Chruslinska et al. 2019; Santoliquido et al. 2021). In addition, constraining the high-redshift peak to z ≳ 8 would provide support for the existence of Pop III binaries.

We first show the inferred distribution of the low-redshift peak, z_l , for all universes in Figure 4. All measurements of z_l constrain the low-redshift peaks to z < 3, with 95% credible interval uncertainties of ∼40%. The true values are contained within the uncertainty.

Zoom In Zoom Out Reset image size

Next, we look at the measurements of the high-redshift peak's location z_h for the universes with Pop III binaries, as shown by the orange violins in Figure 5. In all cases, the true value ${\hat{z}}_{h}=11.6$ lies within the 95% credible intervals. The lower bound of the credible interval is above z ∼ 7 for all values of ${\tilde{f}}_{{\rm{G}}}$, indicating that one can confidently place the secondary peak at a redshift much larger than where the star formation peaks. The widths of the credible intervals decrease from ∼50% to ∼25% when ${\tilde{f}}_{{\rm{G}}}$ increases from zero to 1. This may again be explained by the steeper redshift tail in the cluster population, which makes the Pop III peak easier to resolve.

Zoom In Zoom Out Reset image size

3.2. Phenomenological Analysis

Having shown how a simple nonparametric model can already provide insight into the existence of a high-redshift population of BBHs, we now repeat the hierarchical inference using the phenomenological model described in Appendix C.

We start by showing the posteriors on the peak of the merger rate density for the Pop III BBHs, this time obtained as one of the parameters of the phenomenological Pop III model, ⁷ in Figure 5 (green violins). Remembering that the orange violins in the same plot report the measurement we obtained with the GPR approach, we find that the two methods yield very consistent results, with the widths of the 95% credible intervals varying by ∼10% at most. The consistency between the two approaches highlights the promise of the nonparametric approach in revealing the existence and location of a high-redshift population. However, the phenomenological model directly describes the morphology of each subpopulation and thus allows for extracting information about each individual population, which cannot be accessed with the nonparametric approach.

In Figure 6, we show the inferred $\dot{n}$ for the simulated ${\tilde{f}}_{{\rm{G}}}=0.5$ universes with (right panel) and without (left panel) Pop III binaries. ⁸ For each population, the solid line represents the true merger rate, whereas the dashed line and colored band represent the median and 68%/95% credible intervals.

Zoom In Zoom Out Reset image size

The total $\dot{n}$ (black band) can be well constrained within a few percent level up to z ∼ 6. For f_III = 0, the relative uncertainty rises to ∼100% at z ≳ 10. Even at low redshift, the uncertainty of ${\dot{n}}_{{\rm{F}}}$ (orange band) and ${\dot{n}}_{{\rm{G}}}$ (blue band) is about 50%, ∼10 times larger than that of the total rate, $\dot{n}$. This is because the morphology of ${\dot{n}}_{{\rm{F}}}$ and ${\dot{n}}_{{\rm{G}}}$ is similar at z ≲ 3, where most of the BBHs can be detected with a precise distance measurement. Conversely, ${\dot{n}}_{{\rm{F}}}$ and ${\dot{n}}_{{\rm{G}}}$ are easier to distinguish at z ≳ 5, where, as discussed before, ${\dot{n}}_{{\rm{G}}}$ is declining more rapidly than ${\dot{n}}_{{\rm{F}}}$, which instead has a long tail at z ≳ 8. Overall, the similarity in the morphology of the low-redshift volumetric merger rate of the two dominating channels induces degeneracies, and hence boosts uncertainties, in the individual merger rates, ${\dot{n}}_{{\rm{F}}}$ and ${\dot{n}}_{{\rm{G}}}$.

Considering now the universe with f_III = 0.024 (right panel), we find that rate of the Pop III curve (green band) near its peak at z ∼ 12 can be measured with a relative uncertainty of ∼50%, while the uncertainty of other channels remains similar to the universe with f_III = 0.

We compare the phenomenological recovery of the total $\dot{n}$ to the GPR recovery (Figure 2) for the same universes and find that the typical uncertainty of $\dot{n}$ is smaller by a factor of ∼2. This is not surprising, since the phenomenological approach uses models for the subpopulations, which inform the recovery of the overall merger rate. Naturally, the price to pay for the improved precision is to have made the results depend on the goodness of the models.

It is worth looking at the correlations between the hyperparameters of the various subpopulations. Some correlation should be expected, since, for example, the number of sources at high redshift might be potentially explained either by the model with a larger fraction of field binaries, which have a fat high-redshift tail, or by binaries in the Pop III channel. In Figure 7, we show the marginalized 2D contours of the $({\tilde{f}}_{{\rm{G}}},{f}_{\mathrm{III}})$ pair for the universes ${\tilde{f}}_{{\rm{G}}}=0.5$ with (yellow) and without (purple) Pop III. A positive correlation between the two parameters is clearly visible. This is caused by the partial model degeneracy between ${\dot{n}}_{{\rm{F}}}$ and ${\dot{n}}_{\mathrm{III}}$. This goes exactly in the direction one would expect; underestimating f_III means that the model must increase the number of field binaries to account at least partially for the high-redshift binaries. But if the number of field binaries increases, ${\tilde{f}}_{{\rm{G}}}$ must decrease.

Zoom In Zoom Out Reset image size

The partial correlation between ${\dot{n}}_{{\rm{F}}}$ and ${\dot{n}}_{\mathrm{III}}$ manifests itself in two other interesting ways. First, we observe an increase in the uncertainty of ${\tilde{f}}_{{\rm{G}}}$ when Pop III mergers are present. In Figure 8, we show violin plots for the marginalized ${\tilde{f}}_{{\rm{G}}}$ posteriors at different true ${\tilde{f}}_{{\rm{G}}}$ values. The yellow and purple violins correspond to the universes with and without Pop III binaries, respectively. While the true values lie inside the 95% credible interval in all cases, the uncertainties increase by ∼10% for the universes with Pop III binaries.

Zoom In Zoom Out Reset image size

Second, the uncertainty in f_III decreases as ${\tilde{f}}_{{\rm{G}}}$ increases. Figure 9 is a plot similar to Figure 8 but showing the marginalized f_III posteriors in different universes. When Pop III binaries are present, the uncertainty stays roughly constant, Δf_III ∼ 0.1, up to ${\tilde{f}}_{{\rm{G}}}\lesssim 0.5$, after which it drops gradually to Δf_III ∼ 0.05 at ${\tilde{f}}_{{\rm{G}}}=1$. A similar trend is observed when no Pop III binaries exist; the uncertainty drops from ∼2.5% to ≲0.1%. This is because ${\dot{n}}_{{\rm{F}}}$ has a longer tail in the high redshift z ≳ 8 and is easier to confuse with a small excess of Pop III.

Zoom In Zoom Out Reset image size

In this paper, we have shown that observations made by a network of 3G detectors can be used to infer the properties of different BBH populations using redshift-only information. The larger horizon of a 3G detector network allows for accessing thousands of BBHs per month up to z ∼ 15, which is necessary to resolve the excess of high-redshift (z ≳ 8) BBHs originating from Pop III stars. We consider ∼16,000 binaries, roughly corresponding to 2 months of data, and multiple values of the branching ratio between binary formation in galactic fields and globular clusters. For every value of this branching ratio, we analyzed both a case where a few hundred Pop III BBHs are present and one where they do not exist.

First, we consider a hierarchical inference approach based on a nonparametric reconstruction of the total volumetric merger rate density $\dot{n}$. We look for local peaks in the reconstructed total merger rate and extract limited but useful information about the high-redshift population. By requiring that a possible high-redshift peak have an amplitude of at least 1/10 of the low-redshift peak, we find evidence for the presence of a high-redshift peak when Pop III binaries are included and constrain its position to be in the range 7 ≤ z_h ≤ 15 for various mixing fractions between field and cluster binaries. Using the same approach, we rule out the existence of secondary peak structure if there were no Pop III binaries. This minimally modeled measurement of the position of a high-redshift peak (or lack thereof) in the total binary merger rate might, with some models, be translated into a measurement or an upper limit on the abundance of Pop III stars. With a similar approach, we are able to measure the position of a low-redshift peak, which might be used to investigate typical time delays between start formation and mergers.

Then, we considered a modeled analysis where a phenomenological model exists for each of the three subpopulations, which are characterized by a set of unknown hyperparameters measured from the data together with the (unknown) branching ratios. Among the most remarkable results, we found that, irrespective of the true value of the relative abundances of field and cluster binaries, the Pop III fraction can be constrained to be f_III ≳ 0.01 (f_III ≲ 0.02) at 95% credibility for the universes that have (do not have) Pop III binaries. In both cases, the branching ratio between field and cluster binaries can be measured with better than ∼100% uncertainty (95% credible interval). The precision of the measurement of Pop III mainly depends on the morphology of the merger rate densities of the dominating channels in the high-redshift region. If the dominating channels predict a shallower declining slope at high redshift, an eventual contribution to the high-redshift merger rate from Pop III is less distinctive, introducing correlations with the dominating channels.

Some studies suggest that Pop III might contribute to a nonnegligible fraction of the merger rate in the local universe (Hartwig et al. 2016; Liu & Bromm 2020a) or even a secondary peak in the low redshift z ∼ 2 due to different formation scenarios (Kinugawa et al. 2020; Liu & Bromm 2020b). If this additional low-redshift peak exists, it will make the Pop III subpopulation more distinguishable while degrading the measurements of branching ratios, owing to extra degeneracies in the low-redshift regime.

We emphasize that our analysis is only assuming two months worth of data. A back-of-the envelope calculation assuming that the statistical uncertainty shrinks like $1/\sqrt{N}$ would imply a factor of ∼5 improvement over the results we present here after 5 yr of data. In that scenario, the phenomenological approach could identify a fraction of Pop III mergers as small as f_III ∼ 0.5%. However, as stressed multiple times in this work, the phenomenological inference requires reliable models for the merger rate of the various formation channels and will yield results that are as good as the models. On the other hand, the model-independent approach, though intrinsically less precise, has the attractive feature of not requiring any specific modeling of the underlying subpopulation. In this analysis, we used a three-detector network. A smaller network would lead to worse redshift measurements for individual sources and hence yield worse statistical uncertainties than reported in this paper. However, this can be compensated for by a longer observation time.

In this work, we have only considered three subpopulations: the galactic field and cluster binaries and high-redshift Pop III binaries. As mentioned above, many other channels have been proposed and can plausibly contribute a sizable fraction of the total rate. Our analysis can be trivially extended to include these and other subpopulations at the price of increasing computational cost and correlations and potentially degrading the measurement of some of the parameters. On the other hand, most of these different channels predict distinctive features in the BBHs they produce besides their redshift distribution, for example, masses, spins, and eccentricity (Dominik et al. 2012, 2013, 2015; de Mink & Belczynski 2015; Rodriguez et al. 2015, 2016; Belczynski et al. 2016; Bartos et al. 2017; O'Shaughnessy et al. 2017; Raidal et al. 2017, 2019; Stevenson et al. 2017; Vitale et al. 2017; Rodriguez & Loeb 2018; Di Carlo et al. 2019; Mapelli et al. 2019; Yang et al. 2019; Yi & Cheng 2019; Antonini & Gieles 2020; Biscoveanu et al. 2020; Breivik et al. 2020; Gröbner et al. 2020; Kremer et al. 2020; Samsing et al. 2020; Santoliquido et al. 2021; Yang et al. 2020; Tagawa et al. 2020a, 2020b, 2021). Including these features can enhance the precision of multipopulation inference, help break degeneracy between population models, and improve the understanding of each formation channel. We leave this extension of multidimensional BBH parameters in the 3G era to a future work.

The authors would like to thank Emanuele Berti, Hsin-Yu Chen, Carl Haster, and the LVK's rate and population working group for fruitful discussions and comments. The authors also thank Katarina Martinovic, Carole Perigois and Tania Regimbau for cross-checking the validity of the Pop III phenomenological model. K.N. and S.V. acknowledge the support of the National Science Foundation through the NSF award PHY-1836814. K.N. and S.V. are members of the LIGO Laboratory. W.M.F. is funded by the Center for Computational Astrophysics at the Flatiron Institute, which is supported by the Simons Foundation. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-1764464. This paper carries LIGO document No. LIGO-P2000540 and CE document No. CE-P2000007.

Here we briefly review the main statistical tool used in the analysis, i.e., the hierarchical Bayesian inference framework (Farr et al. 2015; Mandel et al. 2019; Thrane & Talbot 2019; Wysocki et al. 2019; Vitale 2020). We model the production mechanism of BBHs as an inhomogeneous Poisson process whose differential merger rate in the detector frame ⁹ is given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}^{d}}{{dz}}(z| {\boldsymbol{\Lambda }},{R}^{d})=\displaystyle \frac{1}{1+z}\displaystyle \frac{{dR}}{{dz}}(z| {\boldsymbol{\Lambda }},R),\end{eqnarray} \tag{ A1 }$$

where dR/dz is the differential merger rate in the comoving frame characterized by the "shape parameters" Λ and an overall normalization factor given by the total merger rate in the comoving frame $R=\int \displaystyle \frac{{dR}}{{dz}}{dz}$. The factor 1/(1 + z) accounts for the cosmological time dilation effect so that the total merger rate in the detector frame is

$$\begin{eqnarray}&&{R}^{d}=\int \displaystyle \frac{{{dR}}^{d}}{{dz}}{dz}=\int \displaystyle \frac{1}{1+z}\displaystyle \frac{{dR}}{{dz}}{dz}.\end{eqnarray} \tag{ A2 }$$

The vector of shape parameters Λ contains the quantities that are used to model the underlying physical populations (see Section D).

One can write the hyperposterior of the population parameters Λ and R^d given a set of N_obs observations ${\boldsymbol{d}}\equiv {\{{d}_{i}\}}_{i=1}^{{N}_{\mathrm{obs}}}$ as

$$\begin{eqnarray}&&\begin{array}{l}p\left({\boldsymbol{\Lambda }},{R}^{d}| {\boldsymbol{d}}\,\right)\\ \quad \propto \left[\prod _{i=1}^{{N}_{\mathrm{obs}}}{T}^{d}\displaystyle \int {{dz}}_{i}\,p\left({d}_{i}| {z}_{i}\right)\displaystyle \frac{{{dR}}^{d}}{{dz}}\left({z}_{i}| {\boldsymbol{\Lambda }},{R}^{d}\right)\right]{e}^{-{R}^{d}{T}^{d}}\\ \,\times \,\pi \left({\boldsymbol{\Lambda }},{R}^{d}\right)\\ \quad \simeq \left[\prod _{i=1}^{{N}_{\mathrm{obs}}}{T}^{d}\displaystyle \frac{1}{{M}_{i}}\sum _{j=1}^{{M}_{i}}\displaystyle \frac{{{dR}}^{d}}{{dz}}\left({z}_{{ij}}| {\boldsymbol{\Lambda }},{R}^{d}\right)\right]{e}^{-{R}^{d}{T}^{d}}\\ \qquad \times \,\pi \left({\boldsymbol{\Lambda }},{R}^{d}\right).\end{array}\end{eqnarray} \tag{ A3 }$$

In going from the second to the last line of Equation (A3), we have approximated the integrals with discrete sums. For the ith source, this amounts to calculating an average of the merger rate evaluated at the M_i points ${\{{z}_{{ij}}\}}_{j=1}^{{M}_{i}}$ drawn from the likelihood p(d_i ∣z_i ) of the ith source. In the third line, π(Λ, R^d ) is the hyperprior, and T^d is the experiment duration in the detector frame.

We find it more useful to quote the volumetric merger rate density, defined for the kth subpopulation as

$$\begin{eqnarray*}&&{\dot{n}}_{k}(z)\equiv {{dR}}_{k}/{{dV}}_{c},\end{eqnarray*}$$

where dV_c /dz is the differential comoving volume. Then, the merger rate history ${{dR}}_{k}^{d}/{dz}$ in the detector frame as a function of redshift for the kth subpopulation is

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}_{k}^{d}}{{dz}}(z| {{\boldsymbol{\Lambda }}}_{k},{R}_{k}^{d})=\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{\dot{n}}_{k}(z| {{\boldsymbol{\Lambda }}}_{k},{R}_{k}),\end{eqnarray} \tag{ A4 }$$

where the subscript k denotes the relevant quantities of the kth subpopulation.

The overall merger rate can be expressed as the sum of the individual merger rates of all P subpopulations (P = 3 in our analysis), i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}^{d}}{{dz}}(z| {\boldsymbol{\Lambda }},{{\boldsymbol{R}}}^{d})=\sum _{k}^{P}\displaystyle \frac{{{dR}}_{k}^{d}}{{dz}}(z| {{\boldsymbol{\Lambda }}}_{k},{R}_{k}^{d}),\end{eqnarray} \tag{ A5 }$$

where the vectors Λ and R ^d contain all Λ_k and ${R}_{k}^{d}$, respectively. Since ${\sum }_{k}^{P}{R}_{k}=R$ and ${\sum }_{k}^{P}{R}_{k}^{d}={R}^{d}$, we may rewrite Equation (A5) in terms of the branching ratios in the detector frame, i.e., the fraction of sources in each subpopulation, ${f}_{k}\equiv {R}_{k}^{d}/{R}^{d}$,

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}^{d}}{{dz}}(z| {\boldsymbol{\Lambda }},{R}^{d})={R}^{d}\sum _{k}^{P}{f}_{k}{p}_{k}^{d}(z| {{\boldsymbol{\Lambda }}}_{k}),\end{eqnarray} \tag{ A6 }$$

where ${p}_{k}^{d}(z| {{\boldsymbol{\Lambda }}}_{k})$ is the normalized merger rate of the kth population in the detector frame, and the f_k are subject to the constraint ${\sum }_{k}^{P}{f}_{k}=1$. Since we expect the fraction of Pop III binaries to be small, it is more convenient to introduce the fraction of cluster binaries over the sum of field and cluster binaries, ${\tilde{f}}_{{\rm{G}}}\equiv {f}_{{\rm{G}}}/({f}_{{\rm{G}}}+{f}_{{\rm{F}}})$.

Therefore, for the parameterized analysis, we model the merger rates of the three subpopulations in terms of several phenomenological parameters, which are treated as unknowns, together with the branching ratios (Appendix D). The total merger rate in Equation (A3) is thus calculated by adding up the contribution of each channel. On the other hand, in the unmodeled approach, we measure the overall $\displaystyle \frac{{dR}}{{dz}}\,(z| {\boldsymbol{\Lambda }},R)$ directly, without making any assumption about the individual subpopulations that might be contributing to it.

More details on the two approaches can be found in Appendices B and C, whereas the hyperpriors are described in Appendix E. Throughout the study, we assume Planck 15 cosmology (Ade et al. 2016).

This section provides details on the implementation of the GRP that we use to infer $\dot{n}$ without assuming any specific functional form. We only require that dR/dz is sufficiently smooth such that dR/dz can be described by a piecewise-constant function over W = 30 redshift bins, which are uniformly distributed in linear space in the range 0 ≤ z ≤ 15. The merger rate dR/dz is thus written as

$$\begin{eqnarray}\displaystyle \frac{{dR}}{{dz}}=\left\{\begin{array}{ll}\displaystyle \frac{{\rm{\Delta }}{R}_{1}}{{\rm{\Delta }}{z}_{1}} & 0\leqslant z\lt {z}_{1}\\ \ldots & \\ \displaystyle \frac{{\rm{\Delta }}{R}_{i}}{{\rm{\Delta }}{z}_{i}} & {z}_{i-1}\leqslant z\lt {z}_{i}\\ \ldots & \\ \displaystyle \frac{{\rm{\Delta }}{R}_{W}}{{\rm{\Delta }}{z}_{W}} & {z}_{W-1}\leqslant z\lt {z}_{W}\end{array}\right.,\end{eqnarray} \tag{ B1 }$$

where ΔR_i is the merger rate in the ith redshift bin Δz_i ≡ z_i − z_i−1 = 0.5 so that ${\sum }_{i=1}^{W}({\rm{\Delta }}{R}_{i})\,\equiv \,R$. To make the GPR more efficient, we infer dR/dz in natural-log space. Then, we apply a squared-exponential Gaussian process prior on ${X}_{i}\equiv \mathrm{ln}({\rm{\Delta }}{R}_{i})$, with a covariance kernel

$$\begin{eqnarray}\begin{array}{rcl}{K}_{{ij}} & \equiv & \mathrm{Cov}\left({X}_{i},{X}_{j}\right)\\ & = & {\sigma }_{X}^{2}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{z}_{i-1/2}-{z}_{j-1/2}}{l}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ B2 }$$

where ${z}_{i-1/2}=\tfrac{1}{2}\left({z}_{i}-{z}_{i-1}\right)$ is the midpoint of the ith redshift bin, ${\sigma }_{X}^{2}$ is the variance of {X_i }, and l is the correlation length in redshift space. The multivariate Gaussian process prior on the random variable vector ${\boldsymbol{X}}\equiv \{\mathrm{ln}({\rm{\Delta }}{R}_{i})\}$ with a mean vector μ _X and covariance matrix K ≡ {K_ij } is then

$$\begin{eqnarray}&&{ \mathcal G }({\boldsymbol{X}}| {{\boldsymbol{\mu }}}_{X},{\sigma }_{X},l)\equiv { \mathcal N }\left[{\boldsymbol{X}}| {{\boldsymbol{\mu }}}_{X},{\boldsymbol{K}}\left({\sigma }_{X},l\right)\right].\end{eqnarray} \tag{ B3 }$$

The kernel K enforces the smoothness of dR/dz on scales that are comparable to or larger than l, which may be much larger than the bin spacing if the data support it, and prevents overfitting when W is large (Foreman-Mackey et al. 2014). To further enhance the sampling efficiency, we utilize the Cholesky factorization to decompose K into a lower-triangular matrix L , such that

$$\begin{eqnarray}&&{\boldsymbol{X}}\equiv {\mu }_{X}+{\boldsymbol{L}}({\sigma }_{X},l){\boldsymbol{\eta }}\,\end{eqnarray} \tag{ B4 }$$

follows the same Gaussian process prior ${ \mathcal G }$ by drawing η from a multivariate standard normal distribution. A common choice of μ _X is a constant mean vector μ_X 1. Since we know that dR/dz has a strong dependence of the redshifted differential comoving volume dV_c /dz/(1 + z), we further impose the mean of the Gaussian process prior to be the natural log of dV_c /dz/(1 + z) (normalized to N_obs within the comoving volume V_c and the observation time T) with a common shift Δμ_X , i.e.,

$$\begin{eqnarray}&&\{{\mu }_{X,i}\}=\left\{\mathrm{ln}\left(\displaystyle \frac{{N}_{\mathrm{obs}}}{{V}_{c}T}{\int }_{{z}_{i-1}}^{{z}_{i}}\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{dz}\right)+{\rm{\Delta }}{\mu }_{X}\right\}\,,\end{eqnarray} \tag{ B5 }$$

and we treat Δμ_X , which is a single variable, as an additional parameter to include any possible fluctuation in the normalization R. We then obtain $\dot{n}$ from dR/dz divided by the differential volume in each bin. All together, the W + 3 = 33 hyperparameters in the GPR are thus

$$\begin{eqnarray}&&{{\boldsymbol{\Lambda }}}_{\mathrm{GPR}}=\left({\{{\eta }_{i}\}}_{i=1}^{W},{\rm{\Delta }}{\mu }_{X},{\sigma }_{X},l\right).\end{eqnarray} \tag{ B6 }$$

Directly modeling the volumetric merger rate $\dot{n}$ of a subpopulation as a function of astrophysical quantities, such as various distributions of initial stellar mass, the mass/radius of star clusters, BH natal kicks, or stellar metallicity, requires detailed stellar evolution or N-body simulations that are computationally expensive. To facilitate our analysis, we model the three formation channels phenomenologically. For the field, we follow the Madau–Dickinson functional form, ¹⁰

$$\begin{eqnarray}&&{\dot{n}}_{{\rm{F}}}(z| {\alpha }_{{\rm{F}}},{\beta }_{{\rm{F}}},{C}_{{\rm{F}}})\propto \displaystyle \frac{{\left(1+z\right)}^{{\alpha }_{{\rm{F}}}}}{1+{\left(\tfrac{1+z}{{C}_{{\rm{F}}}}\right)}^{{\beta }_{{\rm{F}}}}},\end{eqnarray} \tag{ C1 }$$

where α_F, β_F, and C_F are unknown parameters that characterize the upward slope at z ≲ C_F − 1, the downward slope at z ≳ C_F − 1, and the peak location of the volumetric merger rate density, respectively.

For cluster binaries, we describe the volumetric merger rate as a lognormal distribution in cosmic time, which we treat as a function of redshift,

$$\begin{eqnarray}&&{\dot{n}}_{{\rm{G}}}(z| {\mu }_{{\rm{G}}},{\sigma }_{{\rm{G}}},{t}_{{\rm{G}}})\propto \mathrm{LogNorm}(t(z)-{t}_{{\rm{G}}}| {\mu }_{{\rm{G}}},{\sigma }_{{\rm{G}}}),\end{eqnarray} \tag{ C2 }$$

where t(z) is the cosmic time as a function of redshift, and LogNorm is the standard lognormal distribution of the argument t − t_G parameterized by μ_G and σ_G. The additional parameter, t_G, is a reference time that marks the birth of the first cluster binaries.

For Pop III, we use the following functional form:

$$\begin{eqnarray}&&{\dot{n}}_{\mathrm{III}}(z| {a}_{\mathrm{III}},{b}_{\mathrm{III}},{z}_{\mathrm{III}})\propto \frac{{e}^{{a}_{\mathrm{III}}(z-{z}_{\mathrm{III}})}}{{b}_{\mathrm{III}}+{a}_{\mathrm{III}}{e}^{({a}_{\mathrm{III}}+{b}_{\mathrm{III}})(z-{z}_{\mathrm{III}})}},\end{eqnarray} \tag{ C3 }$$

where a_III, b_III, and z_III characterize the upward slope at z < z_III, the downward slope at z > z_III, and the peak location of the volumetric merger rate density, respectively.

We have verified that these three phenomenological models can well fit the data from population synthesis analysis (Belczynski et al. 2016, 2017; Rodriguez & Loeb 2018) for the values of their arguments given in Equations (C1)–(C3).

We define the branching ratio between field and cluster binaries, ${\tilde{f}}_{{\rm{G}}}$, implicitly through the equations

$$\begin{eqnarray}&&{f}_{{\rm{G}}}\equiv {\tilde{f}}_{{\rm{G}}}\left(1-{f}_{\mathrm{III}}\right),\end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray}&&{f}_{{\rm{F}}}\equiv \left(1-{\tilde{f}}_{{\rm{G}}}\right)\left(1-{f}_{\mathrm{III}}\right),\end{eqnarray} \tag{ C5 }$$

where {f_k } are the original fractions of Equation (A6).

Therefore, there are a total of 12 hyperparameters in the phenomenological model:

$$\begin{eqnarray}&&{{\boldsymbol{\Lambda }}}_{\mathrm{PM}}=\left({\alpha }_{{\rm{F}}},{\beta }_{{\rm{F}}},{C}_{{\rm{F}}},{\mu }_{{\rm{G}}},{\sigma }_{{\rm{G}}},{t}_{{\rm{G}}},{a}_{\mathrm{III}},{b}_{\mathrm{III}},{z}_{\mathrm{III}},{\tilde{f}}_{{\rm{G}}},{f}_{\mathrm{III}},{R}^{d}\right).\end{eqnarray} \tag{ C6 }$$

In this section, we describe how we prepare the simulated universes that will be analyzed with the methods described in the previous section.

First, we need to choose reference (i.e., "true") merger rate densities that will be used to generate the redshift of the BBHs. We do so by means of the phenomenological curves in Equations (C1)–(C3). As described in Appendix C, these curves describe the morphology of each volumetric merger rate density and are parameterized by, e.g., the rising slope at low redshift, the declining slope at high redshift, and the redshift at which the merger rate peaks. The phenomenological curves are obtained by fitting Equations (C1)–(C3) to the simulation results available in the literature (Belczynski et al. 2016, 2017; Rodriguez & Loeb 2018). Specifically, we take the model-averaged simulation results of Belczynski et al. (2016) and Rodriguez & Loeb (2018) for field and cluster binaries, respectively, as well as the "FS1" model's result of Belczynski et al. (2017) for Pop III binaries. We stress that our phenomenological fits include the effect of the time delay from binary formation to merger, as well as the impact of stellar metallicity on the binary evolution specified in the population synthesis analyses. In particular, this implies a quite remarkable fact: the same Madau–Dickinson function can be used to fit both the SFR (this is its normal use) and the merger rate density it implies for different values of its parameters.

We use the following numbers as the "true" values of each curve parameter when preparing our sources:

$$\begin{eqnarray*}\begin{array}{rcl}\left({\hat{\alpha }}_{{\rm{F}}},{\hat{\beta }}_{{\rm{F}}},{\hat{C}}_{{\rm{F}}}\right) & = & (2.57,5.83,3.36),\\ \left({\hat{\mu }}_{{\rm{G}}},{\hat{\sigma }}_{{\rm{G}}},{\hat{t}}_{{\rm{G}}}\right) & = & (1.63,0.96,0.66\,\mathrm{Gyr}),\\ \left({\hat{a}}_{\mathrm{III}},{\hat{b}}_{\mathrm{III}},{\hat{z}}_{\mathrm{III}}\right) & = & (0.66,0.3,11.6).\end{array}\end{eqnarray*}$$

While this fixes the true shape of the merger rate for each subpopulation, we still need to fix the amplitudes. We use nine different values of the relative merger rate between cluster and field binaries, ${\tilde{f}}_{{\rm{G}}}$, equally spaced in the range from zero to 1. For each value of ${\tilde{f}}_{{\rm{G}}}$, we consider two universes: one with and one without Pop III binaries.

Following the latest GWTC-2 result, we fix the local volumetric merger rate density to $\dot{n}(0)=25\,{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$ (Abbott et al. 2020b). This yields ${R}_{{\rm{F}}}^{d}+{R}_{{\rm{G}}}^{d}\approx 8000$ month^–1. In the universes with Pop III binaries, the Pop III fraction f_III ≈ 0.024 is chosen to generate an additional 200 sources month^–1 coming from the Pop III channel. This number is chosen such that the peak merger rate density of Pop III binaries is ∼10 times smaller than the peak of the dominating channels. The nominal value of this peak, ${\dot{n}}_{\mathrm{III}}({z}_{\mathrm{III}})\approx 20\,{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$, is consistent with the comparison of the "FS1" model to the Population I/II field binaries in Belczynski et al. (2016). We stress that the current predictions for the Pop III merger rate span a few orders of magnitude in either direction relative to the one we are using (Kinugawa et al. 2014, 2016; Belczynski et al. 2016; Hartwig et al. 2016). This is because the formation efficiency of Pop III binaries greatly depends on the initial mass function of the Pop III stars and the distribution of the initial orbital separation (see Belczynski et al. 2016; Hartwig et al. 2016). Since ${\dot{n}}_{\mathrm{III}}$ declines rapidly at later cosmic time, it does not significantly contribute to $\dot{n}(0)$.

To ensure a reasonable measurement of luminosity distance, which typically requires three or more detectors, we choose a baseline detector network with one CE in Australia, one CE in the United States, and one ET in Europe (Vitale 2016; Vitale & Evans 2017; Vitale & Whittle 2018). Generally speaking, not all BBHs within the detector horizon can be detected with S/Ns over some thresholds, depending on their orientation and intrinsic parameters. This results in a Malmquist bias (Mandel et al. 2019; Vitale 2020). However, the efficiency only drops to 50% at z ∼ 20 for a typical 30 M_⊙ and 30 M_⊙ BBH mass-pair. Therefore, we limit our analysis to the redshift range 0 ≤ z ≤ 15 and neglect selection effect.

Finally, we simulate a month of data in the following way. We first draw the set of true redshifts, {z_true}, from the "true" dR/dz in each universe. Then, for each z_true,i , we obtain the observed redshift, z_obs,i , by drawing a random variable from a mock-up single-event likelihood conditional on z_true,i . Following Vitale et al. (2019), we approximate the likelihood for redshift as a lognormal distribution conditional on z_true,i with a standard deviation ${\sigma }_{\mathrm{LN},i}=0.017{z}_{\mathrm{true},i}$. Finally, we draw 100 single-event likelihood samples conditional on z_obs,i with the same ${\sigma }_{\mathrm{LN},i}$ as calculated previously. The second step is necessary in order to generate a scattering to the true value within the probable range of the single-event likelihood function. Otherwise, the alignment of the true value and the mean of the likelihood introduces systematic bias in the analysis.

To summarize, we generate 18 simulated universes. In all universes, we generate 16,000 (2 months worth of data) BBHs from field and cluster binaries with a given branching ratio ${\tilde{f}}_{{\rm{G}}}$ and add 400 more Pop III observations for the nine universes containing Pop III binaries.

The priors of all of the Λ_GPR hyperparameters are tabulated in Table 1.

Table 1. Hyperpriors for the GPR Model

ΛGPR,i	Prior Function	Prior Parameters	Domain
ηi	Normal	(μN, σN) = (0, 1)	(−∞ , +∞)
ΔμX	Normal	(μN, σN) = (0, 10)	(−∞ , +∞)
σX	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=(0,4)$	(0, +∞)
l	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=\left(0,\tfrac{1}{2}\mathrm{ln}(10)\right)$	(0, +∞)

Download table as:  ASCIITypeset image

The priors of all of the Λ_PM hyperparameters are tabulated in Table 2. The choice of lognormal priors is made to restrict each model to always carrying a characteristic peak within the redshift range, rather than increasing or decreasing monotonically. The Gaussian prior and the domain of z_III ensure that the peak of ${\dot{n}}_{\mathrm{III}}$ lies at high redshift z ≥ 8 and prevents ${\dot{n}}_{\mathrm{III}}$ from mimicking the two dominating channels that have peaks at low redshift z ≲ 3.

Table 2. Hyperpriors for the Phenomenological Models

ΛPM,i	Prior Function	Prior Parameters	Domain
αF	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{\alpha }}_{{\rm{F}}},0.25)$ a	(0, 10]
βF	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{\beta }}_{{\rm{F}}},0.25)$	(0, 20]
CF	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{C}}_{{\rm{F}}},0.25)$	(0, 6]
μG	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{\mu }}_{{\rm{G}}},0.25)$	(0, 5]
σG	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{\sigma }}_{{\rm{G}}},0.25)$	(0, 5]
tG	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{t}}_{{\rm{G}}},0.25)$	(0, 2.0]
aIII	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{a}}_{\mathrm{III}},0.5)$	(0, 2]
bIII	Lognormal	$({\mu }_{\mathrm{LN}},{\sigma }_{\mathrm{LN}})=({\hat{b}}_{\mathrm{III}},1)$	(0, 2]
zIII	Normal	$({\mu }_{{\rm{N}}},{\sigma }_{{\rm{N}}})=({\hat{z}}_{\mathrm{III}},2)$ b	[8, 20]
fIII	Uniform	⋯	[0, 0.5]
${\tilde{f}}_{{\rm{G}}}$	Uniform	⋯	[0, 1.0]
Rd	Half Cauchy	γC = Nobs c	[0, +∞]

^a Here ${\mu }_{\mathrm{LN}}$ and ${\sigma }_{\mathrm{LN}}$ are the mean and standard deviation of the lognormal distribution, respectively. ^b Here μ_N and σ_N are the mean and standard deviation of the normal distribution, respectively. ^c Here γ_C is the scale parameter of the Cauchy distribution.

Download table as:  ASCIITypeset image

In this section, we report the population hyperposteriors for the modeled analysis, i.e., the posteriors of the variables that parameterize the individual subpopulations described in Appendix C, as well as the branching ratios.

First, we show the hyperposterior of the dominating channels' parameters for the universe with $({\tilde{f}}_{{\rm{G}}},{f}_{\mathrm{III}})=(0.5,0)$ (the same as Figure 6(a)) in Figure 10. Most of the shape parameters are measured with ∼100% uncertainty, e.g., ${\tilde{f}}_{{\rm{G}}}={0.59}_{-0.27}^{+0.28}$. We note that there are a few interesting correlated pairs, such as $({t}_{{\rm{G}}},{\tilde{f}}_{{\rm{G}}})$ and (t_G, σ_G). Since t_G characterizes the starting time of the merging cluster binaries, an earlier t_G shifts the peak of the cluster population toward a higher redshift, where the cosmological volume is smaller. Hence, ${\tilde{f}}_{{\rm{G}}}$ needs to be larger to keep the same number of cluster binaries. On the other hand, a smaller t_G tends to shift the cluster population toward a lower redshift. Then, a larger σ_G is necessary to maintain a wide merger rate peak.

Zoom In Zoom Out Reset image size

Next, we look at the hyperposterior of the Pop III merger rate parameters for the universe with $({\tilde{f}}_{{\rm{G}}},{f}_{\mathrm{III}})=(0.5,0)$ (the same as Figure 6), as shown by the purple contours in Figure 11. We are able to constrain f_III ≲ 0.02. Both a_III and b_III are very close to their priors, which make sense, since f_III is small, which implies that no information about ${\dot{n}}_{\mathrm{III}}$ can be gained. On the other hand, the position of the peak, z_III, is shifted toward lower values relative to its prior to suppress the contribution from ${\dot{n}}_{\mathrm{III}}$ to the high-redshift total merger rate, which in this universe is entirely determined by the field binaries.

Zoom In Zoom Out Reset image size

For comparison, we also report the results for the universe with f_III = 0.024 (this is the same as in Figure 6(b)). In Figure 11, the yellow contours show the hyperposterior of the ${\dot{n}}_{\mathrm{III}}$ parameters. While the shape parameters of a_III and b_III are still not well constrained, the peak ${z}_{\mathrm{III}}={11.75}_{-1.91}^{+1.92}$ is measured with ∼30% relative uncertainty. Importantly, f_III = 0, i.e., the absence of a Pop III channel, is excluded from the 95% credible interval of the marginalized f_III posterior. This provides strong evidence of the existence of Pop III binaries in our simulated data.