The $\chi_\mathrm{eff}-z$ correlation of field binary black hole mergers and how 3G gravitational-wave detectors can constrain it

Understanding the origin of merging binary black holes is currently one of the most pressing quests in astrophysics. We show that if isolated binary evolution dominates the formation mechanism of merging binary black holes, one should expect a correlation between the effective spin parameter, $\chi_\mathrm{eff}$, and the redshift of the merger, $z$, of binary black holes. This correlation comes from tidal spin-up systems preferentially forming and merging at higher redshifts due to the combination of weaker orbital expansion from low metallicity stars given their reduced wind mass loss rate, delayed expansion and have smaller maximal radii during the supergiant phase compared to stars at higher metallicity. As a result, these tightly bound systems merge with short inspiral times. Given our fiducial model of isolated binary evolution, we show that the origin of a $\chi_\mathrm{eff}-z$ correlation in the detectable LIGO--Virgo binary black hole population is different from the intrinsic population, which will become accessible only in the future by third-generation gravitational-wave detectors such as Einstein Telescope and Cosmic Explorer. Finally, we compare our model predictions with population predictions based on the current catalog of binary black hole mergers and find that current data favor a positive correlation of $\chi_\mathrm{eff}-z$ as predicted by our model of isolated binary evolution.


Introduction
The detection of gravitational waves (GWs) from coalescing binary black holes (BBHs) by the LIGO-Virgo-KAGRA (LVK) collaboration has opened a new window for the study of stellar and binary astrophysics (Aasi et al. 2015;Acernese et al. 2015;Akutsu et al. 2021). To date, the LVK collaboration has reported 69 BBH events with a false alarm rate (FAR) smaller than 1 yr −1 (Abbott et al. 2019(Abbott et al. , 2021a. However, after more than half a decade since the first detection of GWs, the origin of merging BBHs remains an open question. This is not due to a lack of theoretical predictions but rather because of the degeneracy between different formation channel model predictions and unconstrained astrophysical processes of these models (see, e.g., Mandel & Broekgaarden 2022;Zevin et al. 2021).
Improved sensitivity of the LVK detectors and planned thirdgeneration (3G) GW detectors such as the Einstein Telescope (Punturo et al. 2010) and the Cosmic Explorer (Reitze et al. 2019) will increase BBH detection rates by orders of magnitude. A larger sample size allows for detailed investigations of correlations between BBH observable properties (e.g., Maggiore et al. 2020;Tiwari 2021), which might enable different astrophysical formation channels to be distinguished. For example, multiple studies have looked for potential correlations between masses and redshifts (Fishbach et al. 2021;Abbott et al. 2021d), M chirp − χ eff (Safarzadeh et al. 2020;Abbott et al. 2021d;Franciolini & Pani 2022), χ eff − q (Callister et al. 2021b;Abbott et al. 2021d), and χ eff − z (Biscoveanu et al. 2022). The redshift at which the BBH systems merge, z, is a proxy for the distance to the source, M chirp = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 is the chirp mass where m 1 and m 2 are the BH component masses, q = m 2 /m 1 is the mass ratio defined with m 2 < m 1 , and χ eff = (m 1 a 1 + m 2 a 2 )/(m 1 + m 2 ) ·L is the effective spin parameter where a 1 and a 2 are the component BH dimensionless spin vectors andL the orbital angular momentum unit vector.
Here, we demonstrate that field-formed BBHs naturally predict a χ eff − z correlation. Under the assumptions of efficient angular momentum transport inside stars, supported by asteroseismology observational constraints (Kurtz et al. 2014;Deheuvels et al. 2014;Gehan et al. 2018) and current GW observations (Belczynski et al. 2020;Zevin et al. 2021), and Eddington lim-Article number, page 1 of 14 arXiv:2204.02619v2 [astro-ph.HE] 27 Jun 2022 A&A proofs: manuscript no. aanda ited mass accretion efficiency onto BHs, the origin of BH spin in field BBHs arises from tidal interactions during the late BH-Wolf-Rayet (BH-WR) (Qin et al. 2018;Bavera et al. 2020;Fuller & Lu 2022) or WR-WR (Hotokezaka & Piran 2017;Olejak & Belczynski 2021) evolutionary phases or, alternatively, through chemically homogeneous evolution induced by rotational mixing caused from tidal spin-up during the early evolutionary stage of close binaries (Mandel & de Mink 2016;Marchant et al. 2016). From first principles, such correlation should exist since the strength of tidal interaction steeply depends on the orbital separation (Zahn 1977;Hut 1981), and the distribution of orbital separation pre core collapse evolves with redshift. The redshift evolution of the orbital separation distribution originates from the metallicity-dependent stellar winds (Nugis & Lamers 2000;Vink et al. 2001) whose intensity increases as a function of metallicity. Stronger wind-mass loss (during the BH-WR binary evolution phase) widens the binary more efficiently, inhibiting or even canceling the effects of tides. Because the mean metallicity of the Universe decreases as a function of redshift (Madau & Dickinson 2014;Madau & Fragos 2017), we empirically expect an increasing fraction of BBH mergers with highly spinning BH components from tidal spin up as a function of redshift. Additionally, low metallicity stars are more compact at zero-agemain-sequence (ZAMS), expand later in their evolution and have smaller maximal radii during the supergiant phase compare to stars at higher metallicity.
In this paper, we discuss the evolving χ eff distribution as a function of redshift for field BBHs from the common envelope (CE), stable mass transfer (SMT), and the chemically homogeneous evolution (CHE) channels. The paper is structured as follows. First, we introduce our fiducial model and describe how we quantify the χ eff − z correlation in Section 2. In Section 3, we present the redshift evolution of the χ eff distribution for the intrinsic and detectable BBH population as predicted by our model of isolated binary evolution. We then compare our model predictions against the LIGO-Virgo catalog of BBHs. In Section 4, we discuss how potential uncertainties in our model might affect the χ eff distribution of field BBHs and how a possible change in the mixing fraction of the different channels predicted from isolated binary evolution might impact our results. All findings are summarised in Section 5.

The binary black-hole population synthesis model
This study uses the isolated binary evolution model presented in Bavera et al. (2022a), calculated using the POSYDON framework (Fragos et al. 2022), which accounts for BBH formation through the CE, SMT, and CHE channels. It was shown that this model (i) leads to BBH observable properties consistent with the events of the second LIGO-Virgo GW transient catalog (GWTC-2) (Zevin et al. 2021), (ii) have BBH merger rate estimates compatible with observational constraints of GWTC-2 and, now GWTC-3, (Bavera et al. 2021a;du Buisson et al. 2020), (iii) the subpopulation of highly spinning BBHs might explain the observed population of luminous LGRBs across the cosmic history of the Universe (Bavera et al. 2022a), and (iv) does not violate current upper limit estimates of the stochastic GW background (Bavera et al. 2022b).
In contrast to most rapid population synthesis studies, our simulations accurately model the late tidal spin-up phase of the second-born BH and CHE due to rotational-induced mixing of tidally spun-up ZAMS binaries. The former is done by follow-ing the evolution of the binaries from ZAMS up to the formation of the BH-WR systems after the second mass transfer phase with the rapid population synthesis code COSMIC (Breivik et al. 2020) and then uses detailed MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019 BH-WR simulations (Bavera et al. 2021a) to accurately model the final tidal spin-up phase of the BH-WR system up to central carbon exhaustion of the WR star, as done in Bavera et al. (2022a). The detailed BH-WR simulations self-consistently model the angular momentum evolution of the WR star, which is determined by the interplay of tides, WR stellar wind mass loss, and the evolution of the WR stellar structure.
Massive stars in short orbital periods (p < 2 days) at ZAMS with nearly equal masses tidally spin up to be highly rotating, which induces rotational mixing and eventually leads to CHE. Because COSMIC cannot model the parameter space leading to CHE as the code cannot accurately follow the back-reaction on the stellar structure and evolution from rotational-induced mixing, CHE is done by matching ZAMS binary conditions to detailed MESA simulations targeting CHE according to du Buisson et al. (2020), as implemented in Bavera et al. (2022a).
Given the availability of the stellar profile at carbon exhaustion from the MESA simulations, in both cases, the core collapse considers disk formation during the collapse of highly spinning stars. Additionally, we account for mass loss through neutrinos, pulsational pair-instability and pair-instability supernovae (PPISNe & PISNe) (Marchant et al. 2019), and orbital changes resulting from anisotropic mass loss and isotropic neutrinos mass loss (Kalogera 1996), as explained in Appendix D of Bavera et al. (2021a). Because we implement the Fryer et al. (2012) delayed collapse mechanism which assigns zero velocity kicks to collapsing stars with carbon-oxygen cores with masses above 11 M , in practice, we find a statistically small number of systems with χ eff < 0. Alternatively, non-negligible kicks would lead to a more considerable fraction of negative χ eff (see, e.g., Rodriguez et al. 2016;Gerosa et al. 2018;Callister et al. 2021a;Stevenson 2022). For a detailed explanation of the main features and the physical assumptions made in this model, we refer the reader to the extensive discussions in Bavera et al. (2022a).

Detection rates
Merger rates are computed by convolving the redshift and metallicity dependent star formation rate as predicted by the Illustris-TNG simulation (Nelson et al. 2015) with the synthetic catalog of merging BBHs obtained by evolving initial ZAMS conditions at different discrete metallicities with POSYDON. Following the notation of Bavera et al. (2020Bavera et al. ( , 2021aBavera et al. ( , 2022a, the BBH detection rate of a GW detector network can be expressed as a Monte Carlo sum over the synthetic population of merging BBHs, i.e., R det = i, j,k w i, j,k (p det ) yr −1 where w i, j,k is the weighted contribution of a binary k forming at redshift z f,i and merging at redshift z m,k ≡ z k . Here the dummy index j indicates the discrete sum over the 30 simulated log-binned metallicity intervals ∆Z j . The synthetic BBH population is distributed across the cosmic history of the Universe in the center of time bins of size ∆t i = 100 Myr with center the formation redshift z f,i . We chose the time bin size to be small enough to ensure the convergence of our results (see Appedinx D of Bavera et al. 2022a for the details of the calculation).
To compute the BBH detection rate of LIGO-Virgo, we account for the detectors' selection effects, p det , given the source redshift, BH masses, and spins. Here, we assume a GW detector network configuration composed by LIGO Hanford, LIGO Livingston, and Virgo at O3 mid-high/late-low sensitivity (Abbott et al. 2018) with a network signal-to-noise ratio (S/N) threshold of 12 as implemented in Bavera et al. (2021a).
We also consider the sensitivity of the future 3G groundbased GW detector Einstein Telescope. To approximate the BBH detection rate of the Einstein Telescope, we account for detector selection effects given the source redshift and BH masses assuming a theorized noise-sensitive curve ET-D (Hild et al. 2011) as implemented by Barrett et al. (2018) in COMPAS (Team COM-PAS et al. 2022). Here, we assume a conservative S/N threshold of 12 for the Einstein Telescope, similar to what Hild et al. (2011) assumed. In practice, we find that this assumption sets the horizon of a BBH with m 1 = m 2 = 15 M at z = 10, namely p ET det (z = 10, m 1 = m 2 = 15 M ) 0. Finally, we will distinguish the intrinsic detection rate, i.e., what a GW detector with infinite sensitivity would observe on Earth, using the notationw i, j,k = w i, j,k (p det = 1) as first introduced in Bavera et al. (2022a).

Relative channel contribution
Once detection rates are defined, we can compute the relative contribution to the intrinsic detection rate from each one of the isolated binary evolution channels (CE, SMT, and CHE) at a given redshift as where z ∈ [0, 10] is discredited in bins, ∆z, taken to have a constant cosmic time width of ∆t = 200 Myr. Similarly, for the detectable population, we define f det channel (z) where we use w i, j,k instead ofw i, j,k in z ∈ [0, 1] for LIGO-Virgo and z ∈ [0, 10] for the Einstein Telescope.

Quantifying the χ eff − z correlation
To quantify the redshift evolution of the χ eff distribution, we define f χ eff >χ 0 (z) to be the fraction of merging BBHs with χ eff above the arbitrary value χ 0 at a given redshift for the modeled intrinsic BBH population. This quantity is calculated as and, similarly, for the detectable populations, we define f det χ eff >χ 0 (z) with the same redshift spacing and bounds as in Eq. (1).

Results
We now investigate the χ eff − z correlation of field-formed BBHs and assert its detectability given current and planned GW observatories. We first look at the intrinsic and detectable χ eff distributions as a function of redshift in our fiducial model, which includes potential contribution from the CE, SMT, and CHE channels described in Section 3.1. We then quantify the intrinsic and detectable χ eff − z correlation by computing the quantities f χ eff >χ 0 (z) and f det χ eff >χ 0 (z) and the relative channel contributions f channel (z) and f det channel (z), in Section 3.2. Finally, we look for evidence of the modeled f χ eff >χ 0 (z) and f det χ eff >χ 0 (z) in LIGO-Virgo GWTC-3 data in Section 3.3.

The χ eff distribution of field BBHs
First, we show the χ eff distribution as a function of discrete redshift bins for the intrinsic and detectable BBH populations in Figure 1. At low redshifts, the intrinsic distribution manifests a peak at χ eff = 0 plus an almost flat distribution up to χ eff 0.5 which then progressively decays. Similarly, the detectable LIGO-Virgo χ eff distribution also exhibit a similar narrow peak at χ eff = 0. However, in contrast to the intrinsic distribution, we observe a second broader peak at around χ eff 0.35 with an elongated tail reaching large χ eff depending on redshift. Both distributions evolve with redshift; with the χ eff distribution in the intrinsic population showing a slow evolution with redshift, while in the LIGO-Virgo observable population the distribution evolves significantly over the redshift range between 0 and 1. The medianχ eff value of the intrinsic distribution grows fromχ z∈[0,1] The origin of the redshift evolution of the χ eff distribution is different between the intrinsic and the detectable LIGO-Virgo BBH populations as they probe different redshift horizons, z ∈ [0, ∞] and z ∈ [0, 1], respectively. The former encapsulates all merging BBHs at any redshifts and probes the increasing fraction of systems experiencing tidal spin up prior to BBH formation at increasing redshifts (see Section 1). In contrast, the detectable LIGO-Virgo population is biased by the BBH massdependent selection effects. More massive BHs can be detected at further distances than lighter BHs. We note that at the highest redshifts detectable by LIGO-Virgo (z 1), only systems with large positive χ eff are detected due to the increased duration of the inspiral and therefore the S/N. However, the impact of χ eff on detectability is minor compared to mass selection effects (Ng et al. 2018).
In the following section, we show how the different evolutionary channels CE, SMT, and CHE, which have distinct spin distributions, have different detector horizons due to their inherent mass spectrum. For a discussion about the intrinsic and LIGO-Virgo observable joint distributions of χ eff vs. M chirp , we refer the reader to Figure 1 of Bavera et al. (2022a). Finally, because the Einstein Telescope has a much more distant horizon than current generation GW detectors, the planned GW observatory will be able to detect the majority of the underlying BBH population up to large redshifts, see e.g., z ∈ [4, 5] in Figure 1. We therefore find that the Einstein Telescope will observe an evolving χ eff distribution similar to the intrinsic one.

The χ eff − z correlation of field BBHs
The χ eff − z correlation of field BBHs in the intrinsic population, f χ eff >χ 0 (z), is shown in the leftmost column of Figure 2 for χ 0 = 0.2 and χ 0 = 0.5. In both cases, f χ eff >χ 0 (z) is monotonically increasing and reaches an asymptotic plateau at high redshifts, z > 5 and z > 8 for χ 0 = 0.2 and χ 0 = 0.5, respectively. To understand the origin of the f χ eff >χ 0 (z) shape, we need to look at this quantity channel-wise and consider the relative contribution of each channel f channel (z) to the total BBH intrinsic population. In Figure 2, we can see that at low redshifts the intrinsic BBH merging population is composed of a mix of channels, f CE (z = 0) = 30%, f SMT (z = 0) = 55%, and f CHE (z = 0) = 15%. In contrast, at higher redshifts, the total population of merging BBHs is dominated by the CE channel, with f CE (z ≥ 2) ≥ 80%.  Hild et al. (2011). In all cases, the fraction of non-spinning BBHs decreases as a function of redshift, shifting the χ eff distribution to larger χ eff values.

Fig. 2.
Fractions of BBHs with f χ eff >0.2 and f χ eff >0.5 as a function of redshift (solid lines), and the relative contribution of each field BBH channel, f channel , according to the legend (dashed lines). Black solid lines show the fraction of systems that satisfy a given χ eff criteria for the combined CE, SMT, and CHE channels with redshift-dependent branching fractions accounted for. (Left) We show the modeled intrinsic (underlying) population of field merging BBHs. (Center) We show the modeled LIGO-Virgo detectable BBH population assuming simulated O3 detector sensitivity selection effects. (Right) We show the modeled Einstein Telescope detectable BBH population assuming forecast detector sensitivity as in Hild et al. (2011). In most cases, the fraction of highly spinning BBHs increases as a function of redshift.
Channel-wise, we can see that f CE χ eff >χ 0 (z) increases monotonically due to a larger fraction of systems experiencing tidal spinup as a function of redshift. On average, at higher redshifts, binary systems are born at lower metallicities and experience reduced stellar wind mass loss. Hence, an increased fraction of binaries can maintain short orbital separations and tidal lock-ing during the BH-WR phase. A similar argument can be made for the SMT channel. However, because SMT leads on average to wider BH-WR orbital separations than the CE channel (Bavera et al. 2021a); we find f SMT χ eff (z)>χ 0 (z) < f CE χ eff >χ 0 (z) for any redshift. Moreover, we find that at low redshifts f SMT χ eff >0.2 (z 0) = 0, which steadily increases to f SMT χ eff >0.2 (z 10) 0.5. On the contrary, f CHE χ eff >0.2 (z) manifests a monotonically decreasing behavior. At low redshifts, f CHE χ eff >0.2 (z 0) = 1, i.e. all BBHs from the CHE channel are fast spinning, while at higher redshifts most systems possess negligible χ eff , with f CHE χ eff >0.2 (z = 10) = 0. We also notice that highly rotating CHE systems with χ eff > 0.5 are not present in the local universe f CHE χ eff >0.5 (z = 0) = 0 but their presence peaks at f CHE χ eff >0.5 (z = 2.5) = 0.6 before decreasing again to f CHE χ eff >0.5 (z = 10) = 0. The decreasing fraction of highly rotating CHE systems as a function of redshift is a direct consequence of angular momentum loss due PPISNe. The CHE channel only operates at low metallicities (Z < 5 · 10 −3 ) but only binaries with metallicities Z ≤ 10 −4 experience mass loss due to PPISNe (in the considered ZAMS primary masses range ≤ 150 M ). During these pulses, the mass ejection from the stellar surface depletes the angular momentum content of these stars and leads to slowly rotating BHs, see Appendix A for more details. Because at large redshifts (z > 3) the metallicity-dependent star formation rate leads to an increasing relative fraction of extremely low metallicity binaries, we expect to observe a decreasing f CHE χ eff >χ 0 as a function of increasing redshift. The three channels combined lead to the monotonically increasing behavior of f χ eff >χ 0 (z) we see in Figure 2, which is mainly dominated by the CE channel. In Appendix B, we show how our fiducial model predictions would change if one of these three channels would be neglected, see Section 4 for a discussion of these alternative scenarios.
The χ eff − z correlation of field BBHs in the detectable populations, f det χ eff >χ 0 (z), is shown in the center and right columns of Figure 2 for LIGO-Virgo detectors at O3 sensitivity and the Einstein Telescope, respectively, for χ 0 = 0.2 and χ 0 = 0.5. For LIGO-Virgo detectability, f det χ eff >χ 0 (z) is a monotonically increasing function growing from f det χ eff >0.2 (z = 0) = 0.25 to f det χ eff >0.2 (z = 1) = 1, and f det χ eff >0.5 (z = 0) = 0.05 to f det χ eff >0.5 (z = 1) = 0.5. On the other hand, the Einstein Telescope f det χ eff >χ 0 (z) mimics the intrinsic distribution up to z 5 above which it shows a suppression. The similarity between the Einstein Telescope detectable distribution and the underlying distribution is due to the Einstein Telescope GW horizon being much more distant than that of LIGO-Virgo. The suppression for the detectable Einstein Telescope f det χ eff >χ 0 (z > 5) is caused by the fact that the detector cannot resolve all distant low mass highly rotating BBHs formed through the CE channel. To understand the difference between the LIGO-Virgo f det χ eff >χ 0 redshift evolution compared to the Einstein Telescope and the intrinsic BBH population, we need to once again consider the relative contribution of each channel. Similar to the intrinsic distribution, for LIGO-Virgo at low redshift, we have a mixed contribution of the different channels, Up to redshift z = 0.4 the SMT channel dominates over CE and CHE, above which the CHE channel dominates the LIGO-Virgo detectable population to the point where f det CHE (z ≥ 0.75) 1. This is notably different than the behavior of the intrinsic BBH population and is a direct consequence of selection effects favouring high BH masses. The different channels have different BH mass distributions, which result in different observational horizons for each channel. Notably, the CHE channel leads to more massive BBHs compared to CE and SMT, and hence this channel can be probed by LIGO-Virgo at larger redshifts compared to BBHs formed from the CE and SMT channels. The described signature leads to a bimodal distribution of χ eff in the LIGO-Virgo detectable BBH population in Figure 1. .2 and f det χ eff >0.5 are obtained from the median of 10,000 GWTC-3 mock catalog events obtained by sampling the 69 events with FAR < 1 yr −1 likelihoods. The modeled prediction for O3 sensitivity is shown with a solid orange line. To compare the model with the data, we generated 10,000 mock catalogs of 69 events, to which we added mock measurements uncertainties. We indicate the median and 90% CI modeled fractions with orange dashed line and shaded area, respectively. Mock uncertainties are obtained from the zero-centered GWTC-3 event likelihoods.
The second peak is mostly composed of BBHs formed through the CHE channel (see Bavera et al. 2022a, for further discussions). This bimodal feature is not present in the χ eff distribution of the intrinsic BBH population for z ∈ [0, 1] (see the right panel of Figure 1). Because intrinsically f CHE χ eff >0.2 (z ≤ 1) 1, f SMT χ eff >0.2 (z ≤ 1) 0, and f CE χ eff >0.2 (z) is monotonically increasing at low redshifts, we also find a monotonically increasing f det χ eff >0.2 (z) function for LIGO-Virgo sensitivity. A similar argument can be made for f det χ eff >0.5 (z), where for low redshifts, it holds that f SMT χ eff >0.5 (z ≤ 1) 0 while both f CHE χ eff >0.5 (z ≤ 1) and f CE χ eff >0.5 (z ≤ 1) are monotonically increasing, which results in f det χ eff >0.5 (z) monotonically increasing. Finally, notice that the CHE dominance is not present in the Einstein Telescope detectable population as the 3G detector, given our theorised sensitivity curve, will be able to observe the entire intrinsic BBH population up to redshift z 4 − 5.
3.3. Evidence for the χ eff − z correlation in GWTC-3 data Our model provides a falsifiable prediction that both the underlying and detected high-χ eff fractions f χ eff >χ 0 (z) and f det χ eff >χ 0 (z) for the O3 LIGO-Virgo detector network should increase as a function of redshift if isolated evolution channels dominate the BBH merger rate at low redshifts. We notice that at low redshifts, z < 1, the evolution of this fraction for the intrinsic BBH population is mild, but can be amplified by the selection effects of current ground-based detectors. We use BBH events from GWTC-3 to infer both f χ eff >χ 0 (z) and f det χ eff >χ 0 (z) and compare them against our model predictions. As in Abbott et al. (2021d), we only consider GWTC-3 events with a false alarm rate (FAR) smaller than 1 yr −1 . In GWTC-3, there are 76 events satisfying this condition from which we exclude the binary neutron stars (NSs) GW170817 and GW190425_081805, the NS-BH systems GW190426_152155, GW200105_162426, GW200115_042309, and the events GW190814, GW190917_114630 in which the less massive compact objects have masses that could be either a massive NS or a BH. Our BBH sample therefore includes a total of 69 BBH events.
We first approximate the observed f det χ eff >χ 0 (z) for O3, or f GWTC−3 χ eff >χ 0 (z), in a model agnostic way directly from the observed events. We measure f GWTC−3 χ eff >χ 0 (z) on a sample of 10,000 mock GWTC-3 catalogs composed of 69 BBH events. A mock catalog of events, {x i } N=69 i=1 , is generated by drawing a 2D sample, x k i = (χ k eff , z k ) i , from each event's x i 2D posterior distribution 1 p(χ eff , z|x i ) weighted by the inverse of the prior 2D probability density p(χ eff , z) in order to sample from the likelihood. The events' posterior and prior distributions are publicly released by the LIGO-Virgo collaboration. We approximate the discretelysampled prior distribution probability density function (PDF) with a 2D kernel density estimator (KDE) trained on the GWTC-3 event samples where the bandwidth of the KDE is set by Scott's rule (Scott 2015) as implemented in the Gaussian KDE function of the SciPy Python module (Virtanen et al. 2020). The accuracy of our KDE method to represent the inferred 2D distributions is verified by comparing the histogram of the original samples and samples generated from the KDEs.
To perform a fair comparison of our model with the observations, we need to account for (i) the statistical variance of drawing a sample of 69 events from our model and (ii) to account for the measurement uncertainty for BBH parameters. This is done by generating 10,000 mock samples of 69 events from our model, to which we add mock uncertainty to each event. We approximate measurement uncertainties following a procedure first shown in Bavera et al. (2020) for the χ eff parameter. Here, we extend this procedure to the 2D case. Mock uncertainties are obtained by shifting another set of 10,000 mock GWTC-3 catalogs by each event's median valuex i = (χ eff ,z) i . When the mock uncertainty is added to the model mock samples, we find that this methodology overestimates the measurement uncertainty of events with low redshift of merger. This occurs because this methodology does not assign smaller measurement uncertainties to events with smaller redshifts of merger. Such correlation is expected because of the larger measurement uncertainty for more distant events, which is due to their typically smaller S/N 1 In contrast to the GWTC-3 official analysis, for events in O3b we use posterior samples from the IMRPhenomXPHM analysis as the Mixed and SEOBNRv4PHM analyses do not come with associated prior samples in the 8th November 2021 data release (LIGO Scientific Collaboration, Virgo Collaboration and KAGRA Collaboration 2021c). compared to events merging at lower redshifts. In practice, we find that this procedure only leads to 0.6% of the sample having nonphysical values |χ eff | > 1 and 4.2% of systems having nonphysical z < 0, which we map back to |χ eff | = 1 or z = 0. We claim that this bias is small and does not affect our results as we still find that events with largerz have broader χ eff distributions.
In Figure 3, we show the comparison of the median f GWTC−3 χ eff >χ 0 (z) computed on the sample of mock GWTC-3 catalogs with the model prediction. The GWTC-3 quantity is independently measured for each mock catalog in the interval z ∈ [0, 1] by counting the events meeting the χ eff > χ 0 condition in the discrete redshift bin ∆z with constant cosmic time bin of size of ∆t = 1.6 Gyr and then quote the median value at each redshift bin. We then compare our model predictions with the inclusion of mock uncertainties by overlaying the median and 90% CI f det χ eff >χ 0 (z). We conclude that our model cannot be ruled out given the current GWTC-3 sample. Even though our model 90% CI overlaps with the median inferred GWTC-3 value, a closer comparison with the model median indicates that our model slightly overpredicts the fraction of highly spinning BBHs. This could be due, e.g., to an overprediction of the fraction of highly spinning BBHs formed from the CHE channel which dominates over BBHs formed from the CE channel in the LIGO-Virgo detectable population (see Appendix B) or the existence of an additional channel contributing to the detectable BBH population with small BH spins (see e.g., dynamical formation in globular clusters, Zevin et al. 2021). Other model uncertainties are discussed in Section 4.
We next infer the underlying fraction f χ eff >χ 0 (z) by fitting a model for the astrophysical BBH population to the GWTC-3 data. We jointly fit the mass (m 1 , m 2 ), spin χ eff , and redshift z distribution, allowing the χ eff distribution to evolve redshift but for simplicity neglecting possible correlations between other parameters: p pop (m 1 , m 2 , χ eff , z) = p(m 1 , m 2 )p(χ eff | z)p(z). ( For the mass distribution, p(m 1 , m 2 ), we use the Broken Power Law model from Abbott et al. (2021b) and for the redshift distribution, p(z), we assume the merger rate evolves as a power law in (1 + z) (Fishbach et al. 2018). We model the redshiftdependent spin distribution p(χ eff | z) as a mixture model between a "zero-spin" component, approximated as a narrow Gaussian centered at χ eff = 0 with standard deviation 0.03, and a "positive spin" component, for which we use a Gaussian distribution N T with mean 0.2 < µ p < 0.5 and standard deviation 0.05 < σ p < 0.5 truncated to the range [0, 1] to reflect our model predictions. We take the mixture fraction A between the zero and positive spin components to be a logistic function of z (so that it is always within 0 < A < 1), described by two free parameters, A(z = 0) and A(z = 1). We therefore have where High-χ eff fractions in the underlying distribution inferred from the population fit described in Section 3.3, in blue and orange according to the legend. Lighter contour colors indicate larger CIs of 50% and 90%, respectively. The fraction of BBH systems with large positive spins in the underlying population may increase with increasing redshift (credibility 82%), consistent with our model predictions (black). We do not yet have enough BBH events at z ∼ 1 to accurately measure the χ eff distribution at high z and therefore cannot confidently conclude that the distribution is evolving.
samples for the GWTC-3 BBH events (LIGO Scientific Collaboration, Virgo Collaboration and KAGRA Collaboration 2019, 2020, 2021a,b,c). We use flat priors on all parameters within their ranges specified above. The inferred intrinsic χ eff population distribution at two redshifts, z = 0 and z = 1, is shown in Figure 4. At z = 0 the positive-spin component is constrained to be small, whereas at z = 1, the data permit a larger fraction of systems with high χ eff , although the overall constraints are more uncertain and more closely resemble the prior. We can directly compare the underlying high-χ eff fractions f χ eff >0.2 and f χ eff >0.5 inferred under this fit to the low-redshift z < 1 predictions in the leftmost panel of Redshift and effective spin parameters of the 69 confident BBH observations drawn from the GWTC-3 posteriors (orange; "observed") compared to 69 draws from the inferred distribution fit (blue; "predicted") described in Section 3.3. Each marker shape corresponds to a different set of 69 draws. We plot ten total sets. The inferred model sometimes over-predicts the largest observed χ eff , while the bulk of both observed and predicted draws cover an equivalent portion of the z − χ eff plane in a comparable abundance, confirming that the inferred model is a good fit for the data. Figure 2. In Figure 5, we show the inferred intrinsic f χ eff >0.2 (z) and f χ eff >0.5 (z) versus our astrophysical model predictions. The intrinsic fractions are broadly consistent with the model predictions, although the data prefer slightly smaller fractions of large positive χ eff at all z, similar to the conclusions of Figure 3 regarding the observed fractions. To reiterate, this could be due, e.g., to an overprediction of the contribution of the CHE channel (see Appendix B) or the non-negligible contribution of an additional channel with small BH spins. Other model uncertainties are discussed in Section 4.
We verify the goodness-of-fit of the inferred model by performing posterior predictive checks. Figure 6 shows the comparison between the (z, χ eff ) parameters of ten mock GWTC-3 catalogs versus ten sets of 69 events drawn from the inferred model. This test is analogous to the posterior predictive check in Figure  2 of Fishbach et al. (2021). Each of the ten sets (plotted with a different marker size) corresponds to one draw from the inferred population hyperposterior. We reweight the single-event posterior from each GWTC-3 event to the population distribution specified by the hyperposterior draw, and draw one (z, χ eff ) sample per event. We then draw a set of 69 predicted events from the same population distribution, conditioned on detection. We can see that the inferred model sometimes over-predicts the largest observed χ eff , while the bulk of both the observed and predicted draws cover an equivalent portion of the z − χ eff plane in a comparable abundance, confirming that the inferred model is a good fit to the data.
Despite the suggestive hint that f χ eff >0.2 (z) increases with z, we are not yet able to confidently identify that the χ eff distribution varies with redshift under our parameterization. More precisely, we constrain f χ eff >0.2 (z = 0.3) > 0.06 at 99% credibility and find that f χ eff >0.2 (z) increases with increasing redshift at 82% credibility. Our conclusions are consistent with the results of Biscoveanu et al. (2022), who find that the width of the χ eff distribution likely broadens with increasing redshift, but do not find compelling evidence that the mean χ eff evolves with red-shift. 2 Our model for field BBH formation predicts that the mean of the χ eff distribution must increase with redshift.
We further caution that our phenomenological population fit makes the simplifying assumption that the mass distribution is independent of spin and redshift, despite the fact that we predict a correlation between total mass, χ eff , and redshift. However, given current statistical uncertainties on the inferred intrinsic f χ eff >χ 0 (z), we do not expect our systematic errors on this inferred quantity from mismodeling the population distribution to be significant. However, with future data it will be important to allow for χ eff to vary with both mass and redshift in BBH population fits, because as Figures. 1 and 2 show, some of the observed χ eff evolution in the LIGO-Virgo catalog will be due to an underlying correlation between χ eff and mass. These conclusions are corroborated by Biscoveanu et al. (2022), who find that the preference for χ eff to correlate with redshift is stronger than a possible correlation with primary mass, although the two scenarios can be confused for each other.

Discussion
In this work, we considered a fiducial model for isolated binary evolution. However, model uncertainties can potentially alter BBH observable distributions and rates (see e.g. Broekgaarden et al. 2021 for an extended overview of such uncertainties). Here, we are interested in astrophysical uncertainties which may alter the χ eff − z joint distribution.
Our fiducial model assumes efficient angular momentum transport inside stars which leads to the formation of nonspinning first-born BHs for the CE and SMT channels. Alternatively, a less efficient angular momentum transport would lead to non-negligible birth spins (see e.g. some model variations in Belczynski et al. 2020) which would consequently raise our estimated fraction of highly spinning BBHs f χ eff >χ 0 (z) as the CE and SMT channels dominate the intrinsic BBH population. Nevertheless, current observations are consistent with low birth spins of 0.1 for isolated BHs (Abbott et al. 2021b,d;Miller et al. 2020;Zevin et al. 2021).
In Bavera et al. (2021a), the impact of mass-transfer physics uncertainties on the χ eff distribution of BBHs formed from the CE and SMT channels was investigated, accounting for uncertainties in (i) the unknown efficiency of CE ejection in the α CE −λ parametrization (see, e.g., Ivanova et al. 2013, for a review), (ii) the SMT accretion efficiency onto BHs, and (iii) the criteria for mass-transfer stability. The first uncertainty directly impacts the relative fraction of highly rotating BBHs in the CE channel as the α CE parameter approximately linearly scales with the orbital separation post CE. For a wide range of α CE ∈ [0.2, 5], Bavera et al. (2021a) showed that the BBH fraction of systems with χ eff > 0.1 formed from the CE channel can vary from 0.54 to 0.82 where the merger rate density might also vary by up to one order of magnitude. Nevertheless, Bavera et al. (2021a) showed how both α CE extremes include a non-zero fraction of tidally spun-up BBHs in the CE channel. The second uncertainty affects the initially negligible spin of the first-born BH of a BBH systems formed through the SMT channel. In the case of highly super-Eddington accretion efficiency onto BHs, Bavera et al. (2021a) showed how a non-negligible fraction of first-born BHs could be spun up due to accretion. However, in such cases, depending on the super-Eddington accretion efficiency, Bavera et al. (2021a) found a suppression of the SMT merger rate density up to two orders of magnitude. This occurs because conservative mass transfer is less efficient than unconservative mass transfer in leading to tight BH-WR systems, leading to less BBH systems that can merge in a Hubble time. Finally, the third uncertainty directly impacts the relative fraction of systems that undergo either stable or unstable mass transfer and, hence, CE or SMT evolution. We now examine how these uncertainties might affect the presented χ eff − z correlation.
In the present study, we assumed inefficient CE ejection, namely the model with α CE = 0.5 of (Bavera et al. 2021a). A smaller value than what was assumed here would lead to a more significant fraction of tidal spun-up BBHs. Because the CE channel dominates the intrinsic BBH population, such a scenario would increase the predicted quantity f χ eff >χ 0 (z). In contrast, a more efficient assumption for CE ejection would lead to a smaller fraction of systems that are tidally spun up. For the detectable fraction f det χ eff >χ 0 (z), we expect a small impact of this assumption as the detectable population of BBHs is dominated by the SMT and CHE channels. Since Bavera et al. (2021a) showed that at α CE = 5 there is still a fraction of highly rotating BBHs formed from the CE channel with a medianχ CE eff = 0.16, we can claim that our model will always display a monotonically increasing f χ eff >χ 0 (z) regardless of the α CE value in the CE parameterization.
Our fiducial model assumed Eddington limited mass-transfer accretion efficiency onto BHs. A super-Eddington accretion efficiency onto BHs would boost the fraction of highly spinning BBHs formed from the SMT channel, and, hence, positively contribute to larger values of f χ eff >χ 0 (z) and f det χ eff >χ 0 (z). However, given that the BBH merger rate from the SMT channel (both detected and intrinsic) drops by up to two orders of magnitude compared to our fiducial model when increasing the allowed accretion rate onto BHs (see Table 1 of Bavera et al. 2021a), we would expect a smaller intrinsic contribution to the SMT channel than the one modeled here.
Both uncertainties (i) and (iii) might lead to a smaller relative contribution of the CE channel to the total BBH population than what is assumed here, where the CE channel dominates the f χ eff >χ 0 (z) behavior. Moreover, recent studies employing detailed binary simulations point towards an overestimation of systems evolving through and surviving CEs due to envelope stripping during the CE ceasing earlier than what is assumed in rapid population synthesis codes (Fragos et al. 2019;Quast et al. 2019;Klencki et al. 2021;Marchant et al. 2021;Gallegos-Garcia et al. 2021). Therefore, it is natural to ask ourselves what would happen to the modeled f χ eff >χ 0 (z) fraction if the CE channel is negligible compared to SMT and CHE. In such a scenario, given our model, one would expect that most binaries evolving through CE would either evolve through SMT or successfully emerge from the CE at wider orbital separations. In the first case this would lead to a SMT contribution that is similar or greater than what is modeled here. The second case would lead to a reduced fraction of tidally spun-up CE systems, similarly to the outcome of choosing an efficient α CE values. If the remaining SMT and CHE channels retain a similar fraction of highly spinning BBHs to what is modeled here, one would find f χ eff >0.2 (z < 4) 0.2 which would eventually decay at larger redshifts while the LIGO-Virgo detectable population would still exhibit a monotonically increasing behaviour since the contribution of CE systems to the LIGO-Virgo detectable population is small ( f det CE (z > 0.25) < 10%). In Appendix B, we show how Fig. 2 would change given the omission of the CE channel from our fiducial model. At low redshifts, z < 1, we find that the intrinsic fraction f χ eff >0.2 of this alternative model is still consistent with the 90% CI of GWTC-3 constraints shown in Figure 5.
Similar to this last point, we alternatively entertain the idea of what would happen to the χ eff − z correlation if either the SMT or CHE channels contributions are negligible. This might happen, for example, if we overestimate the contribution of the initial conditions parameter space at low orbital periods that leads to SMT or CHE evolution. In Appendix B, we show that the presented correlation would still be observed. In both cases we still recover monotonically increasing fractions f χ eff >χ 0 (z) and f det χ eff >χ 0 (z). However, we note that the model with the omission of the SMT channel manifests a larger, relatively constant f χ eff >0.2 (z < 1) 0.45 which is inconsistent with the 90% CI of GWTC-3 constraints in Figure 5 of f χ eff >0.2 (z < 0.4) < 0.3. Finally, we find that the model that excludes the CHE channel has both an intrinsic f χ eff >χ 0 and LIGO-Virgo detectable f det

Conclusions
In this paper, we investigated the χ eff − z correlation of fieldformed merging BBHs. An increasing fraction of highly spinning BBHs as a function of redshift is expected. At higher redshifts, stars are formed at lower metallicities, experience weaker stellar wind mass loss, and consequently can maintain their short orbital separations and undergo tidal spin up. We quantified this correlation by the fraction of systems with χ eff > χ 0 as a function of redshift, f χ eff >χ 0 (z). For our fiducial model of field BBHs, which includes the potential contribution of CE, SMT, and CHE channels, this quantity for χ 0 ∈ [0.2, 0.5] shows a monotonically increasing behavior as a function of redshift in the underlying BBH population. We also presented predictions for the detectable f det χ eff >χ 0 (z) for the LIGO-Virgo detector network at O3 sensitivity and the Einstein Telescope. Because of the smaller horizons of current GW detectors (z 1), the origin of the monotonically increasing LIGO-Virgo f det χ eff >χ 0 (z) quantity is different than the intrinsic BBH population or that which the Einstein Telescope will observe in the future. Such differences originate from different BH mass distributions of the various channels. On average, highly rotating BBHs formed from the CHE channel are more massive than tidally spun up BBH systems formed from the CE channel. Hence, LIGO-Virgo detector selection effects favour high BH masses and lead to different observational horizons for different channels. We find that, in contrast to the intrinsic distribution where the χ eff − z correlation is dominated by tidal spun-up BBHs from the CE channel, the CHE channel dominates the LIGO-Virgo detected χ eff − z correlation above z > 0.4.
Finally, assuming isolated binary evolution dominates the detected population of merging BBHs, we performed a model comparison between our fiducial model and LIGO-Virgo GWTC-3 data. We find that current observations favor the prediction of our model that there is a positive correlation between χ eff and z. Such a conclusion is consistent with the results of Biscoveanu et al. (2022) who found that the width of the χ eff distribution likely broadens with increasing redshift, event though they did not find compelling evidence in favor of a redshift evolving mean χ eff . Additionally, our model prediction at low redshifts of a large zero-spin BBH population with an additional subpopulation of systems with spin vectors preferentially aligned to the orbital angular momentum is in agreement with Roulet et al. (2021) and Galaudage et al. (2021) reanalysis of GWTC-2 events. Moreover, our results are consistent with the findings that investigated field BBH observable properties and rates (Bavera et al. 2020(Bavera et al. , 2021a, multi-channel model selection with GWTC-2 data (Zevin et al. 2021), potential constraints from LGRBs (Bavera et al. 2022a), and the current upper limits of the stochastic GW background (Bavera et al. 2022b).
Considering future 3G GW detector facilities, we demonstrated that if isolated binary evolution plays a dominant role in the formation of merging BBHs in the Universe, 3G GW detectors will observe more of the merging BBHs in the Universe and a χ eff − z correlation that is more indicative of the behavior of the underlying population.

Appendix A: Angular momentum loss due to pulsational pair-instability supernovae
Mass loss due to PPISNe can play a role in depleting the angularmomentum reservoir of a collapsing star. Because the pulsations carry away the outer layers of the stars that carry most of the angular-momentum content of the star, this phenomena could have a major impact in reducing the spins of massive BHs.
The impact of PPISNe on the spin of the second-born BH of tidally spun up BH-WR systems was briefly discussed in Zevin & Bavera (2022). For tidally spun-up systems with orbital periods p < 1 day and WR stellar masses of M WR > 40 M at carbon depletion, the first panel of Figure 1 in Bavera et al. (2021b) shows a small suppression of the second-born BH spin obtained from the WR stellar profile collapse of MESA BH-WR simulations from Bavera et al. (2021a). Because WR stellar wind rates scale as a function of metallicity (Vink et al. 2001), only binaries born at low metallicities (prevalently formed at high redshifts) will evolve to have WR stars in such a mass regime. Hence, for the CE channel, we expect this phenomena to have a small impact as on average the channel operates at smaller WR stellar mass. For the SMT channel, we find that in practice this phenomena is relevant only at large redshifts as this channel leads on average to more massive BH-WR star systems compared to the CE channel, resulting in a f SMT χ eff >0.2 (z ≥ 7) 0.45 plateau in Figure 2.
In contrast, we find that the impact of PPISNe onto the spin of BHs formed from the CHE channel is not negligible as this channel only operates at low metallicities (Z < 5 · 10 −3 ) and for massive stars. For metallicities Z ≤ 10 −4 the entire sample of merging BBHs evolving through the CHE channel is formed by stars with ZAMS primary masses 40 M M 1 70 M which undergo PPISN. This occurs because at these low metallicities stellar wind mass loss is weaker compared to larger metallicities, and the stars reach the mass regime of PPISN, see Figure A1 of du Buisson et al. (2020). We note that in our fiducial model we do not simulate BBH formation above ZAMS primary masses of 150 M , hence Figure A1 of du Buisson et al. (2020) should be read accordingly. On the other hand, the 10 −4 < Z ≤ 5 · 10 −3 parameter space leading to the formation of merging BBHs allows for direct collapse and, hence, conservation of angular momentum during the stellar profile collapse (with the exception of extremely highly rotating stars inducing disk formation). In Figure A.1, we show the ZAMS binary conditions leading to merging BBH formation through the CHE channel, showing their final primary BH spins as a function of ZAMS initial orbital period and primary mass which can be directly compared to Figure A1 of du Buisson et al. (2020). We can see that for Z ≤ 10 −4 and ZAMS primary masses 70 M the entire population of BBHs is composed of BBH systems with negligible spins as they have lost their high stellar angular momentum due to PPISN mass ejection. The gap in the parameter space at 1.8 log 10 (M 1 /M ) 2.1 for Z ≤ 10 −4 binaries in Figure A.1 is due to pair-instability supernovae leaving no remnant. For binaries with Z > 10 −4 , this portion of the parameter space is present at larger ZAMS primary masses and orbital periods (see Figure A1 of du Buisson et al. 2020). The impact of PPISN onto the BH spin of BBHs formed from the CHE channel at extremely low metallicities explains the monotonically decreasing behaviour of f CHE χ eff >χ 0 (z) as a function of increasing redshift as the Universe forms more stars at these low metallicities. Distribution of ZAMS binary orbital period, p, primary mass, M 1 , and the final primary BH spin of systems evolving thorough the CHE channel to become merging BBHs in our fiducial model. In this sample we only include BBH systems with inspiral times less than the age of the Universe. Different markers differentiate metallicity regimes according to the legend. For visualisation purposes, we capped the color bar at a BH1 = 0.7 even though there are BHs approaching the general relativistic limit a BH1 = 1. Though binaries with p < 1 day do tidally spin up and evolve through CHE, they later undergo mass loss due to PPISN which depletes the WR star of its angular momentum reservoir.

Appendix B: The χ eff − z correlation with channel exclusion
In this appendix section we show the impact to our results presented in Figure 2 in the hypothetical scenario that one of the three channels considered has a negligible contribution to the formation of merging BBHs.
First, let us consider neglecting the CE channel. Factors that might lead to this hypothetical scenario are discussed in Section 4. Figure B.1 shows how the results presented in Figure 2 would change under this assumption. In this alternative model, the intrinsic fraction f χ eff >0.2 (z) is mainly supported by highly spinning BBHs formed from the CHE channel at z < 5, while at larger redshift the SMT channel contributes with a larger fraction of tidally spun-up BHs. However, we notice that in contrast to our fiducial model the intrinsic fraction f χ eff >0.2 (z) is monotonically decreasing. On the other hand the LIGO-Virgo detectable BBH population shows a similar behaviour as the fiducial model. This occurs as the CE channel contribution to the LIGO-Virgo detectable population is small compared to the SMT and CHE channels, since the CE channel leads to less massive BBHs (cf. Figure 2).
Second, let us consider neglecting the SMT channel. Figure B.2 shows how the results presented in Figure 2 would change under this assumption. Because at low redshifts (z < 5) the SMT channel mostly contributes to the intrinsic distribution with non-spinning BBHs, this alternative scenario leads to a larger f χ eff >0.2 (z) fraction compared to the fiducial model. This hypothetical scenario would result in a LIGO-Virgo detectable Article number, page 11 of 14 A&A proofs: manuscript no. aanda BBH population f det χ eff >0.2 (z) 0.6, in tension with GWTC-3 observations.
Last, let us consider neglecting the CHE channel. As discussed in Section 4 this might occur, for example, in the hypothetical case where the abundance of binary stars at ZAMS with short orbital periods (p < 2 days) is overestimated. This alternative model is presented in Figure B.3. We can see that the f χ eff >χ 0 (z) distribution is similar to what is presented in Figure 2. This is explained by the fact that for any redshift the CHE channel has a small contribution to the intrinsic population of merging BBHs at f CHE (z) < 0.2. On The other hand the LIGO-Virgo detectable f det χ eff >χ 0 (z) manifests an almost flat behaviour up to z = 0.6 above which it sharply increases to reach unity at z 1. This sharp monotonic increase of f det χ eff >χ 0 (z > 0.6) is due to the contribution of tidally spun up BBHs formed from the CE channel completely dominates over BBHs formed from the SMT channel at z > 0.75, as f det CE (z > 0.75) f det SMT (z > 0.75). A comparison between the intrinsic f χ eff >0.2 (z) when excluding one of the three field channels and the inferred distribution given the phenomenological model presented in Eq. (4) is shown in Figure B.4. We can see that a model without the CHE channel is closer to the median inferred intrinsic fraction of f χ eff >0.2 than the fiducial model. Additionally, the model excluding the SMT channel is incompatible with the 90% CI of the inferred fraction as it overpredicts the fraction of highly rotating BBHs.
Article number, page 12 of 14 Bavera et al.: The χ eff − z correlation of field binary black hole mergers and how 3G gravitational-wave detectors can constrain it