The Redshift Evolution of the Binary Black Hole Mass Distribution from Dense Star Clusters

Gravitational-wave detectors are unveiling a population of binary black hole (BBH) mergers out to redshifts z ≈ 1, and are starting to constrain how the BBH population evolves with redshift. We present predictions for the redshift evolution of the BBH mass and spin distributions for systems originating from dense star clusters. Utilizing a grid of 144 state-of-the-art dynamical models for globular clusters, we demonstrate that BBH merger rates peak at higher redshifts for larger black hole primary masses M 1. Specifically, for M 1 ≳ 40 M ⊙, the BBH merger rate reaches its peak at redshift z ≈ 2.1, while for M 1 ≲ 20 M ⊙, the peak occurs at z ≈ 1.1, assuming that the cluster formation rate peaks at z = 2.2. The average BBH primary mass also increases from ∼10 M ⊙ at z = 0 to ∼30 M ⊙ at z = 10. We show that ∼20% BBHs contain massive remnants from next-generation mergers, with this fraction increasing (decreasing) for larger (smaller) primary masses. This difference is not large enough to significantly alter the effective spins of the BBH population originating from globular clusters, and we find that their effective spin distribution does not evolve across cosmic time. These findings can be used to distinguish BBHs from dense star clusters by future gravitational-wave observations.


INTRODUCTION
Gravitational wave (GW) detections of binary black hole (BBH) mergers by Advanced LIGO (LIGO Scientific Collaboration et al. 2015), advanced Virgo (Acernese et al. 2015), and KAGRA (Akutsu et al. 2021) are rapidly expanding our understanding of the formation and evolution of this population of stellar remnants.To date, 90 GW events have been detected by the LIGO-Virgo-KAGRA collaboration (LVK), including ∼ 80 BBH mergers (Abbott et al. 2023a).These detections have sparked searches for their origins, which can offer profound insights into not only the evolutionary processes of massive stars and binaries leading to BBHs but also shed light on their birth environments including dense star clusters (Romero-Shaw et al. 2021;Fishbach & Fragione 2023;Kou et al. 2024).
Limited by statistical uncertainties (due to the small number of detections) and systematical uncertainties (due to imperfect simulations), however, the contributions from different channels are still uncertain, particularly as a function of redshift (e.g., Zevin et al. 2021;Mapelli et al. 2022;Arca Sedda et al. 2023;Cheng et al. 2023;Raikman et al. 2024).A greater number of detections will impose more stringent constraints on the redshift evolution of properties such as the merger rate (e.g., Fishbach et al. 2018), the mass distribution (e.g., Fishbach et al. 2021), and the spin distribution (e.g., Biscoveanu et al. 2022).Furthermore, next-generation GW detectors, such as Einstein Telescope (Punturo et al. 2010), and Cosmic Explorer (Abbott et al. 2017), will be able to probe BBH populations past redshifts z ≳ 10 (Hall & Evans 2019), providing deeper understanding of the formation mechanisms and evolution environments of BBHs.
Here, we study the evolution of the BBH mass distribution as a function of redshift for BBH mergers from globular clusters (GCs), taking into account cluster for-mation properties and rates.Previous studies have only briefly discussed how the BH masses vary at selected redshifts (e.g., Mapelli et al. 2022;Kou et al. 2024;Torniamenti et al. 2024), employed semi-analytical methods to calculate the redshift evolution of merger rates across different BBH chirp masses (Fujii et al. 2017), or focused mostly on the merger rate distributions for high-mass BHs (Rodriguez et al. 2016a).We improve upon earlier studies by utilizing a large grid of GC simulations computed with the state-of-the-art Cluster Monte Carlo code (CMC; e.g., Rodriguez et al. 2022, and references therein) that incorporates various relevant physics for cluster dynamics and the evolution of the stars and compact objects within.We also investigate the detailed characteristics of the redshift evolution of merger rates and masses across BBH mass ranges.
The paper is organized as follows.We describe our methods for calculating the BBH merger rates as a function of redshift from GCs in Section 2. In Section 3, we present the different merger rate distributions corresponding to various BH mass ranges and the BH mass evolution as a function of redshift.We also explore the effective spin distributions as a function of redshift.We discuss the uncertainties of the cluster simulations in Section 4. Finally, we summarize and conclude in Section 5.

Globular Cluster Simulations
We use a large catalog of 144 GC simulations (Kremer et al. 2020) ran with CMC to explore the mass distribution as a function of redshift for merging BBHs from GCs.These catalog models span a wide range of initial conditions, with initial number of stars N = 2 × 10 5 , 4 × 10 5 , 8 × 10 5 , and 1.6 × 10 6 (corresponding to initial masses of about 1.2×10 5 , 2.4×10 5 , 4.8×10 5 , and 9.7×10 5 M ⊙ , respectively), initial virial radius r v = 0.5, 1, 2, and 4 pc, metallicity Z = 0.0002, 0.002, and 0.02, and galactocentric distance r g = 2, 8, and 20 kpc.The galactocentric distance sets the tidal boundaries of the clusters in their host galaxies.All simulations adopt an initial mass function following the standard Kroupa broken powerlaws (Kroupa 2001) in the mass range 0.08 − 150 M ⊙ .The simulations use a King profile (King 1966) with a concentration parameter W 0 = 5 for the stellar distribution.The initial binary fraction is 5% for all primary masses, and the secondary masses are sampled from a uniform mass ratio between 0.1 and 1 of the primary masses (Duquennoy & Mayor 1991).Binary orbital periods are drawn from a distribution flat in log-scale (e.g., Duquennoy & Mayor 1991) ranging from near-contact (≥ 5(R 1 + R 2 ), where R 1 and R 2 are the radii of the bi-nary component stars) to the hard/soft boundary, and the binary eccentricities are drawn from a thermal distribution (e.g., Heggie 1975).The present-day properties of these models match well with those of the observed Milky Way GCs (Kremer et al. 2020, their Figure 2).
The simulations adopt the "rapid model" for stellar remnant formation from core-collapse supernova (Fryer et al. 2012).Specifically, BHs form with mass fallback.Their natal kicks are sampled from a Maxwellian distribution with velocity dispersion of σ = 265 km s −1 (Hobbs et al. 2005), but with reduced magnitude by a factor of 1 − f fb , where f fb is a parameter measuring the mass fraction of the stellar envelope that falls back upon core collapse.
Pulsational-pair instabilities and pair-instability supernovae are also incorporated into the simulations, following the treatment in Belczynski et al. (2016).Massive stars with pre-explosion helium core masses between 45 and 65 M ⊙ undergo pulsations that eject large amounts of mass, leading to a final stellar mass of at most 45 M ⊙ .The collapse of the final products produce BHs of 40.5 M ⊙ (assuming 10% mass loss converting from baryonic to gravitational matter).Stars with helium core masses above 65 M ⊙ are destroyed by pairinstability supernovae leaving behind no remnants.
We assume that all BHs from stellar collapse are born with zero spin.New spins and recoil kicks are applied to the remnants immediately after the BBH mergers, allowing self-consistent computations of BH retention by the model clusters (Rodriguez et al. 2018b, and references therein).CMC includes post-Newtonian dynamics for BHs (Antognini et al. 2014;Amaro-Seoane & Chen 2016;Rodriguez et al. 2018a,b) and models GW emission and BBH mergers during strong encounters selfconsistently.
Figure 1 shows the merger time distributions of various primary mass ranges for all BBHs that merged within a Hubble time from the catalog models.For this plot, we weight all clusters in the CMC catalog equally.Here, t = 0 corresponds to the beginning of the cluster evolution.We find that BBHs with larger primary masses tend to merge at earlier times.In particular, most (∼ 80%) BBHs with M 1 ≳ 40 M ⊙ merged within one Gyr after their host clusters are formed, while only ∼ 30% of BBHs with M 1 ≲ 20 M ⊙ merged in the first Gyr.This is because massive BHs form early in the first few Myrs of stellar evolution.They also mass segregate to the dense cluster cores on similar timescales (Portegies Zwart et al. 2010, their Eq. 20), where frequent gravitational encounters allow them to form binaries and merge at early times.If we approximate the delay time distribution with a power law p(τ ) ∝ τ α , the best-fit power-law indexes are α ≈ −1, −0.8, −1.2, and −1.4 for all primary masses,

Cluster Formation Rate
The population properties of BBHs as a function of redshift depend on their host clusters' initial properties and formation times.We sample GCs according to their initial mass, virial radius, and metallicity.
We assume that the initial cluster mass distribution follows a Schechter function between 10 5 and 10 6 M ⊙ , similar to the present-day masses of Milky Way GCs and corresponding to the mass span of the catalog models.The mass distribution is described by where β S = −2 and the Schechter mass M S = 10 6.26 M ⊙ (Portegies Zwart et al. 2010;Antonini & Gieles 2020a).
For the initial size distribution, we adopt a Gaussian distribution following Fishbach & Fragione (2023) with mean µ r = 2 pc and standard deviation σ r = 2 pc.We truncate the distribution at 0.5 pc and 4 pc1 .We assume a lognormal metallicity distribution at each formation redshift where σ Z = 0.5 and µ Z is given by the following redshift (z) dependent equation from Madau & Fragos (2017), For each redshift, we adopt a cluster formation rate that follows the star formation rate (Madau & Fragos 2017) and is described by where a z = 2.6, b z = 3.6, and z peak = 2.2 is the peak cluster formation redshift.The maximum cluster formation redshift is set to 20, beyond which R GC = 0. Note that while Fishbach & Fragione (2023) has inferred from GW observations that the rate of GC formation may rise more sharply at z > 3 compared to the global star formation rate, simulations of cluster formation indicate that the peak redshift for cluster formation is consistent with being at z ∼ 2 (e.g., Reina-Campos et al. 2022, their Figure 10).Thus, for simplicity, we assume that cluster formation follows the star formation rate.
The cluster formation rate is normalized by the number density of GCs observed at the present day.We define the total number density of dense clusters ever formed as where t max is the maximum lookback time corresponding to a maximum redshift of 20.Many clusters that formed have dissolved or disrupted due to effects including stellar evolution and tidal stripping by their host galaxies.The total number density of dense star clusters formed can be an order of magnitude larger than what we see today (e.g., Antonini & Gieles 2020a).
To take into account the disrupted clusters, we assume n 0 = f disrp n t , where f disrp is the disruption factor and n t is the number density of GCs observed.We adopt n t = 2.31 Mpc −3 (Rodriguez et al. 2016a, and references therein) and adjust the disruption factor such that the total BBH merger rate from GCs at redshift z = 0 is 10 Gpc −3 yr −1 (since the focus of this study is not on the exact amplitude of the BBH merger rate).In this case, f disrp ∼ 5, consistent with Antonini & Gieles (2020a).Note that this disruption factor corresponds to a minimum cluster mass ∼ 10 5 M ⊙ for the catalog models.The BBH merger rate for a cluster is calculated by summing over all BBH mergers from the cluster following the method in Fishbach & Fragione (2023), We weigh each catalog model with its initial mass M i GC , virial radius r j v , and metallicity Z k using the distribution functions and cluster formation rate described in Section 2.2.Here t l is the lookback time, and t m is the time a BBH takes to merge in a simulation.ẑ is the cosmological function converting time to redshift, and ẑ(t l + t m ) is the formation redshift of a cluster.We adopt the cosmological parameters in Planck Collaboration et al. (2016) as implemented in Astropy Collaboration et al. (2018) for the cosmological calculations.The total BBH merger rates from GCs as a function of redshift can be estimated by summing over all catalog models.
Figure 2 shows the BBH merger rates as a function of redshift after taking into account the cluster formation histories.Similar to Figure 1, more massive BHs have their merger rates peak earlier due to mass segregation in GCs and the earlier formation time of massive BHs.In addition, the low metallicity environment at higher redshift also facilitates the formation of more massive BHs (e.g., Torniamenti et al. 2024).The merger rates peak at redshift about 1.1, 1.7, and 2.1 for BBHs with primary masses in the ranges

Star Formation Rate
Figure 2. BBH merger rates from GCs as a function of redshift for four different primary mass ranges.Mergers with larger primary masses have merger rates peaking at earlier times.The peak merger redshifts are 1.2, 1.1, 1.7, and 2.1 with uncertainties of 0.1 for all BBHs, M1 < 20 M⊙, 20 M⊙ ≤ M1 < 40 M⊙, and M1 ≥ 40 M⊙, respectively.The merger rate for all BBHs is normalized to 10 Gpc −3 yr −1 at z = 0 (see also Section 2.2).The cluster formation rate from Madau & Fragos (2017) is also plotted as the dashed curve (arbitrary normalization), with a peak at z = 2.2.
closely the peak of the star and cluster formation rate (z = 2.2) as described by Eq. 5.
The trend of the peak merger redshifts for different BBH primary masses is distinct from that found in BBH mergers from isolated binary evolution.Recently, van Son et al. (2022) showed that common envelope evolution and stable Roche-lobe overflow in isolated binary systems produce BBH mergers with different delay times and masses.BBHs with small primary masses (≲ 30 M ⊙ ) from the common envelope channel dominate the merger rates at high redshift, while BBHs with large primary masses (≳ 30 M ⊙ ) from the stable Rochelobe overflow channel contribute to a large fraction of mergers at low redshift.The combined effects of these two binary evolution processes lead to BBHs with small primary masses merging at higher redshifts than BBHs with high primary masses (van Son et al. 2022, their Figure 8).In particular, they showed that the peak of the merger rates for BBHs with small primary masses M 1 ≲ 10 M ⊙ coincide with the peak of their adopted star formation rate (z = 2.7), while the merger rate for BBHs with M 1 ≳ 30M ⊙ peak at a lower redshift z = 1.9.This is opposite to what we find for BBH mergers produced in GCs as is demonstrated in Figure 2 Figure 3.The evolution of mass as a function of redshift for BBH mergers from GCs.The average primary mass is shown by the dotted orange curve.The dark blue band shows the 25th and the 75th percentiles of the primary mass, while the light blue band shows the 5th and the 95th percentiles.
However, our focus in this work is on the variation of the slopes between different BBH primary masses.We predict the power-law index governing the redshift evolution to vary by ≈ 1 across the mass range.This is too small to be resolvable with the current GW data, which can only constrain the overall merger rate evolution to a 90% credibility width of ≈ 3. Approximately 500 BBH events with the design sensitivity of Advanced LIGO will tighten the 90% credibility width to 1 (Fishbach et al. 2018;Fishbach & Kalogera 2021), so we anticipate that O(1000) events will provide the necessary resolution to identify the predicted variation with BBH primary mass.This will probably require O(1) year of observation with A+ sensitivity (Kiendrebeogo et al. 2023).
We also show the redshift evolution of the primary mass of the BBH population from GCs in Figure 3.The average primary mass increases for larger redshifts, with ⟨M 1 ⟩ reaching about 30 M ⊙ at z ≈ 10 and ⟨M 1 ⟩ ≈ 12 M ⊙ at z ≈ 0. The average mass of BBH mergers at z = 0 is about a factor of two smaller than that estimated by Belczynski et al. (2022), which may be because they used GC models with low metallicities only (stellar winds are weaker at lower metallicities, thus leading to smaller mass loss and more massive BHs; e.g., Belczynski et al. 2010) and they have a smaller number of models.
Recently, Rinaldi et al. (2023) suggested that there is evidence for the evolution of BBH primary mass with redshift from the third GW transient catalog (Abbott et al. 2023a), where the merger rates of BBHs consisting of primary masses ≲ 20 M ⊙ drop at z ≳ 0.4.This may indicate a prevalent contribution of dynamicallyassembled BBHs at higher redshift.However, Abbott et al. (2023b) and Fishbach et al. (2021) showed that within current statistical uncertainties, the GW data is consistent with a nonevolving mass distribution, which matches our expectations that the predicted trend from GCs is not yet resolvable, but will be with future GW observations.

Next-generation Mergers and Spins
In dense star clusters, BBH merger remnants may be retained if the GW recoil velocity is smaller than the escape velocity of the cluster, creating a new generation of BHs.These second-generation (2G) BHs will continue to participate in close dynamical encounters, potentially forming new BBHs that can merge.Therefore, the masses of the merging BBHs may be affected by previous BBH mergers.We define the first-generation BHs (1G) as those formed from collapsing stars and highergeneration BHs as those produced in one or more previous BBH mergers.
Figure 4 shows the primary and secondary masses for BBHs containing different generations of BHs from the catalog models.1G here means both component BHs in a binary are first-generation, while 2G and 3G indicate that at least one BH component is second-or thirdgeneration.Overall, 2G BBHs are more massive than 1G as they contain BH remnants from previous mergers.For primary BH masses in the ranges M 1 < 20 M ⊙ , 20 M ⊙ ≤ M 1 < 40 M ⊙ , and M 1 ≥ 40 M ⊙ , there are about 88%, 83%, and 70% of binaries that consist of only 1G BHs, respectively.
The fractions of BBHs containing next-generation BHs may affect the population's effective inspiral spins as measured with GWs, where ⃗ χ 1 and ⃗ χ 2 are the spin vectors of the BH components, q is the ratio between the secondary and the pri- All 1G Containing 2G Containing 3G mary BH masses, and L is the unit binary angular momentum vector.The 2G BHs created by previous mergers of nonspinning BHs have high spins concentrated at ∼ 0.7 (e.g., Rodriguez et al. 2019), contributing to large effective spins.
To further explore if the BBH effective spins also evolve with redshift, we plot the merger time distributions of different BBH generations and the number ratio of next-generation mergers to 1G mergers in Figure 5.The merger time distribution of 2G BBHs is similar to that of 1G BBHs, except there are no 2G mergers at the very beginning (≲ 20 Myr) of the cluster evolution.The merger timescales of 2G BBHs with M 1 ≳ 40M ⊙ are shorter compared to those of the 1G BBHs with M 1 ≲ 40M ⊙ , whether inside the clusters or following ejection.Thus, despite the slight delay in the merger times of 2G systems compared to their 1G counterparts, the more massive BBHs still merge at higher redshift, as is demonstrated in Figure 2. Overall, the number ratio between 1G and higher-generation mergers stays mostly constant at around 0.2 throughout a Hubble time (consistent with e.g., Rodriguez et al. 2019) and may only start to decline at ≳ 5 Gyr after many BHs have been dynamically processed and ejected out of the clusters.
This indicates that the effective spin distribution of the BBH population from GCs remain mostly unchanged across redshift (assuming BHs are born with zero spins from the collapse of stars).This is distinct from the trend predicted by isolated binary evolution, where it is shown that BBH effective spins tend to increase with increasing redshift due to the tidal spin-up of BBH progenitors at high redshift (e.g., Bavera et al. 2022).Biscoveanu et al. (2022) showed that observationally, the width of the effective spin distribution may broaden with increasing redshift, but found no evolution of the mean effective spin with redshift.Thus BBHs from isolated binary evolution or GCs alone cannot explain this observed trend if BHs are born with zero spin from collapsing stars.In addition, the effective spins also depend on the mass ratios of the component BHs (Eq.8).Most BBH mergers from GCs have mass ratios close to unity regardless of their primary masses (> 50% have q ≥ 0.7 for all M 1 ranges) since mass segregation in dense star clusters leads to the pair-up of similar mass objects.Small mass ratios (q ≲ 0.4; ∼ 10% of all model BBHs) are mostly contributed by BBHs with primary mass M 1 > 20 M ⊙ , and the most extreme mass ratios (q ≲ 0.2; only ∼ 1% of all model BBHs) are produced by BBHs with M 1 > 40 M ⊙ .The limited number of massive BHs, in contrast to the abundance of smaller BHs, often results in the pairing of massive BHs with companions of unequal mass, with the minimum mass ratio constrained by the minimum and maximum BH masses in the models.Nevertheless, the small fraction of BBHs with q ≲ 0.4 in the catalog models makes it improbable for these asymmet-ric systems to significantly affect the effective spins of the overall cluster BBH population.
Only six merging BBHs contain 3G BHs, since 2G BHs created by the previous mergers of two BHs have high spins, and mergers with these fast-spinning 2G BHs produce large GW recoil kicks that will most likely eject the merger remnants from the host clusters (Rodriguez et al. 2019).Note that the very massive 1G BHs in the models (M BH ≳ 100 M ⊙ ) are formed by the collapse of massive stars produced in repeated stellar collisions.

MODEL UNCERTAINTIES AND DISCUSSION
The properties of the BBHs from GCs could be affected by the adopted initial conditions of the clusters and the assumed BH birth spin.
For example, the initial binary fractions of cluster stars are not well constrained.Many GCs at the present day are observed to have low binary fractions ∼ 5 − 10% (Milone et al. 2012).Young massive star clusters in the local universe, on the other hand, are observed to have high binary fractions comparable to the Galactic field (e.g., Sana et al. 2009).In particular, observations have suggested that close to 100% of the O-and B-type stars in the field are born with companions (e.g., Sana et al. 2012;Moe & Di Stefano 2017).Since it is believed that the progenitors of GCs may have similar properties to the young massive star clusters (e.g., Chatterjee et al. 2010Chatterjee et al. , 2013)), massive stars in GCs may have high binary fractions initially.This initial high binary fraction of massive stars can enhance the formation of massive BHs (Arca Sedda et al. 2024) and their mergers (González et al. 2021), in comparison to the catalog models with 5% initial binary fraction for all stars.A larger population of primodial binaries for massive stars may also produce more low-mass BBH mergers (chirp mass ≲ 20 M ⊙ ) at short delay times (Hong et al. 2018), moving the peak redshift of these mergers (e.g., the peak redshift of M 1 < 20 M ⊙ in Figure 2) to higher redshift.However, primordial binaries are unlikely to affect the massive BBH mergers (e.g., Hong et al. 2018, their Figure 9).Furthermore, Chatterjee et al. (2017) showed that primordial binaries containing BH progenitor stars are generally disrupted in GCs due to dynamical encounters, producing BBHs independent of the initial binary properties.Therefore, for dense star clusters like GCs, dynamically-produced BBH mergers probably dominate the merger rates regardless of the initial binary fraction, and the initial binaries may have negligible effects on the mass evolution of BBHs mergers in GCs (Hong et al. 2018;Arca Sedda et al. 2024), preserving the peak redshift trend as a function of the primary BH mass (Figue 2).
The stellar initial mass function can also influence the populations of BBHs originating from GCs. Weatherford et al. (2021) showed that top-heavy initial mass functions (those with higher fractions of massive stars) produce more BBH mergers (Chatterjee et al. 2017), including sources with high BH component masses compared to the canonical initial mass function we adopted (also Wang et al. 2021).
In addition, the birth spins of BHs formed in stellar collapse determine the GW recoil kicks of BBH mergers and the fractions of 2G BHs that can be retained in GCs.Thus, different BH birth spins will also lead to different mass distributions and effective spin distributions of the BBH population from GCs (Rodriguez et al. 2019;Mapelli et al. 2021).
The shapes of the merger rate curves and the absolute merger rates may vary for BBH populations from GCs under different initial conditions and BH birth spin assumptions.However, they are unlikely to alter the trend of the merger rate peaks and slopes as a function of the primary BH mass, as massive BHs will continue to form and undergo mass segregation first.

CONCLUSIONS
We studied the evolution of mass with redshift for BBH mergers produced in GCs.We calculated the BBH merger rates across different primary mass ranges by utilizing a large grid of realistic GC simulations and taking into account both cluster formation rates and properties.
Due to mass segregation in dense star clusters and the early formation of massive BHs, BBHs with massive components tend to merge early in the evolutionary history of their host clusters.The low metallicity environment at high redshift also contributes to the formation of massive BHs.We demonstrated that BBHs with primary masses M 1 ≳ 40 M ⊙ have the peak merger rate at z ≈ 2.1, while BBHs with M 1 ≲ 20 M ⊙ have the peak merger rate at z ≈ 1.1, assuming a cluster formation rate that peaks at z = 2.2.Similarly, the average primary mass increases from ∼ 10 M ⊙ at z = 0 to ∼ 30 M ⊙ at z = 10.Furthermore, we showed that the power-law index κ of the BBH rate evolution at z ≲ 1, presuming R(z) ∼ (1 + z) κ , increases for more massive primary masses from κ ≈ 1 to κ ≈ 2.
The correlation between the redshift of the merger rate peaks and the primary mass of the merging BBHs from GCs is the opposite of that found in isolated binary evolution (van Son et al. 2022).Thus, the redshift evolution of the BBH mass distribution can be used to distinguish their formation channels.Depending on the primary mass spectrum, current GW detections are consistent with a nonevolving mass spectrum but do not yet provide the resolution to probe the degree of mass evolution predicted here (Fishbach et al. 2021).However, we anticipate that during the next couple of observing runs by the LVK, O(1000) BBH detections at z ≲ 2 will provide the constraining power to test our predictions for the mass-dependence of the merger rate evolution.Furthermore, next-generation GW detectors (e.g., Punturo et al. 2010;Abbott et al. 2017) will be able to probe directly past the peak of the merger rate and help investigate the BBH mass distribution across cosmic time.
We also showed that BBH mergers with larger primary masses from GCs contain slightly larger fractions of 2G BHs, assuming BHs are born with zero spins from collapsing stars.About 30% of BBHs with M 1 ≳ 40 M ⊙ consist of at least one 2G component, while only about 12% with M 1 ≲ 20 M ⊙ have 2G components.This is expected since 2G BHs are, in general, more massive than 1G BHs.Despite the increased presence of 2G BHs in the more massive merging BBHs, which also merge at higher redshifts, we observed no redshift-dependent evolution in the effective spins of cluster BBHs.This is due to the overall fraction of BBHs containing 2G BHs remaining relatively constant over a Hubble time.Software: CMC (Joshi et al. 2000(Joshi et al. , 2001;;Fregeau et al. 2003;Fregeau & Rasio 2007;Chatterjee et al. 2010Chatterjee et al. , 2013;;Umbreit et al. 2012;Morscher et al. 2015;Rodriguez et al. 2016bRodriguez et al. , 2022)), Fewbody (Fregeau et al. 2004), COSMIC (Breivik et al. 2020)

Figure 1 .
Figure 1.The distributions of BBH merger times since the formation of their host clusters from the catalog models.Four different primary mass (M1) ranges are shown.About 80% of BBHs with M1 ≥ 40 M⊙ merged within the first Gyr of the evolution of their host clusters, while less than 40% of BBHs with M1 < 20 M⊙ merged in the same period.
and M 1 ≥ 40 M ⊙ , respectively.The merger rate peak of the most massive BHs follows . In the redshift range z ≲ 1 (corresponding to the redshift range of current LVK observations; e.g.Fishbach & van

Figure 4 .
Figure 4.Primary and secondary masses of merging BBHs containing different generations of BHs from the catalog models.There are a total of 9517 BBHs that only contain 1G BHs, 1709 that have at least one 2G BH, and 6 that have at least one 3G BH.

Figure 5 .
Figure 5. Upper panel: the number ratio of next-generation BBH mergers over 1G mergers as a function of merger time since the formation of their host clusters.The ratio stays mostly constant at around 0.2 across a Hubble time except at ≲ 40 Myr when the 2G BBHs have not yet formed.Bottom panel: BBH merger time distributions for BBHs containing three different generations of BHs.
It will likely be larger if lower mass clusters (e.g., ∼ 10 4 M ⊙ ) are taken into account since they produce smaller numbers of BBHs (e.g., Rodriguez & Loeb 2018; Antonini & Gieles 2020b) and dominate the number of clusters formed.