The binary black hole spin distribution likely broadens with redshift

The population-level distributions of the masses, spins, and redshifts of binary black holes (BBHs) observed using gravitational waves can shed light on how these systems form and evolve. Because of the complex astrophysical processes shaping the inferred BBH population, models allowing for correlations among these parameters will be necessary to fully characterize these sources. We hierarchically analyze the BBH population detected by LIGO and Virgo with a model allowing for correlations between the effective aligned spin and the primary mass and redshift. We find that the width of the effective spin distribution grows with redshift at 98.6% credibility. We determine this trend to be robust under the application of several alternative models and additionally verify that such a correlation is unlikely to be spuriously introduced using a simulated population. We discuss the possibility that this correlation could be due to a change in the natal black hole spin distribution with redshift.


INTRODUCTION
The growing catalog of compact binary mergers detected with gravitational waves (Abbott et al. 2019a(Abbott et al. , 2021a has allowed for increasingly precise characterization of the population properties of binary black holes (BBHs) and neutron stars (Abbott et al. 2021c), which can shed light on how these systems form and evolve (Wong et al. 2020;Zevin et al. 2021;Bouffanais et al. 2021). The latest data from the Advanced LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) observatories has revealed that the BBH mass distribution has substructure (Tiwari & Fairhurst 2021;Edelman et al. 2022;Li et al. 2021;Veske et al. 2021) beyond a smooth power-law and single peak or break at ∼ 40 M Talbot & Thrane 2018;Abbott et al. 2019b). BBH spins are found to be small but non-zero (Wysocki et al. 2019;Roulet & Zaldarriaga 2019;Miller et al. 2020;García-Bellido et al. 2021;Biscoveanu et al. 2021), with some of their tilts misaligned to the orbital angular momentum (Talbot & Thrane 2017;Abbott et al. 2019bAbbott et al. , 2021c, although these conclusions have been challenged when applying a different population model (Roulet et al. 2021;Galaudage et al. 2021). The merger rate is found to evolve with redshift at a rate consistent with the star formation rate (Madau & Dickinson 2014;Fishbach et al. 2018;Abbott et al. 2021c).
In addition to fitting the mass, spin, and redshift distributions independently, previous works have looked for correlations among these parameters (Safarzadeh et al. 2020;Callister et al. 2021;Abbott et al. 2021c;Tiwari 2022;Franciolini & Pani 2022). Such correlations can be imprinted via evolutionary processes within a single formation channel or can be caused by the presence of multiple populations arising from distinct formation channels. For example, systems formed dynamically in dense environments via repeated mergers are expected to be more massive and have higher spins (Portegies Zwart & McMillan 2002;Benacquista & Downing 2013;Miller & Lauburg 2009;McKernan et al. 2012;Antonini & Rasio 2016;Rodriguez et al. 2015;Gerosa & Berti 2017;Rodriguez et al. 2019;Kimball et al. 2020;Gerosa & Fishbach 2021). Callister et al. (2021) and Abbott et al. (2021c) found statistically significant arXiv:2204.01578v2 [astro-ph.HE] 9 Aug 2022 evidence for a correlation between the distribution of the mass-weighted spin aligned with the orbital angular momentum, χ eff , and the binary mass ratio, q, where the mean of the χ eff distribution increases for more extreme mass ratios. This sort of correlation is unexpected for most individual BBH formation models, with the possible exception of formation in the disks of active galactic nuclei (McKernan et al. 2012;Stone et al. 2017;Mckernan et al. 2018;McKernan et al. 2020;Tagawa et al. 2020;Callister et al. 2021) and super-Eddington accretion during stable mass transfer for systems formed via isolated binary evolution , so may hint at the superposition of multiple populations formed via independent channels with unique q vs. χ eff signatures. Safarzadeh et al. (2020) examined whether the effective spin distribution correlates with various mass parameters of the binary. They found a possible negative correlation between the mean effective spin and each of the chirp mass (M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 ), total mass, and primary mass parameters and a possible positive correlation between the dispersion of the effective spin and mass; neither finding reached high statistical significance (∼ 80% credibility depending on the mass and correlation parameters). Abbott et al. (2021c) and Tiwari (2022) also found that the spread in the component spins aligned with the orbital angular momentum may increase with the binary chirp mass. This is explained as an effect of the paucity of events at large chirp masses, which manifests as a degradation of the constraint on the spin distribution, rather than a firm measurement of an increase in the width of the distribution at larger chirp masses. Franciolini & Pani (2022) also find initial evidence of an evolution of the χ eff distribution with increasing total mass, which they interpret as evidence for a subpopulation of dynamically-formed or primordial black hole binaries (De Luca et al. 2020a,b,c;, as both models predict correlations between mass and spin. Previous works have also explored potential correlations between the BBH mass and redshift distributions Abbott et al. 2021c). Heavier black holes are predicted to form from lower-metallicity stellar progenitors at higher redshifts, which could lead to such a correlation (Belczynski et al. 2020;Mapelli et al. 2019;Neijssel et al. 2019). For systems formed via isolated binary evolution, those that undergo a common envelope event tend to be both less massive and have shorter delay times, merging at higher redshifts (van Son et al. 2021). Motivated by these theoretical predictions, Fishbach et al. (2021) find that redshift depen-dence of the maximum black hole mass is required if the primary mass distribution has a sharp cutoff, but a gradual tapering of the primary mass distribution is consistent with no evolution with redshift, a finding which was confirmed with the latest catalog of events by Abbott et al. (2021c).
In this work, we search for correlations between the χ eff distribution and the primary mass and redshift distributions. We find robust evidence of a correlation between χ eff and redshift, where the width of the χ eff distribution increases with redshift. We also find a weaker correlation between the width of the χ eff distribution and primary mass. When allowing the χ eff distribution to correlate with both redshift and primary mass, we find a preference in the data for a correlation with one or the other, although we cannot distinguish which. In Section 2, we describe the Bayesian methods and models employed in our analysis. The results on data from the latest GWTC-3 catalog (Abbott et al. 2021b) of compact binaries observed by LIGO-Virgo are presented in Section 3. Section 4 includes a validation of our results using simulated populations and alternative models applied to the real data. We conclude and comment on the potential astrophysical implications of our finding in Section 5.

METHODS
We employ the framework of hierarchical Bayesian inference to constrain the hyper-parameters governing the population-level distributions of the masses, spins, and redshifts of the binary black hole systems detected by LIGO and Virgo. We assume the distribution of primary masses, m 1 , is described by the sum of a truncated power law with low-mass smoothing and a Gaussian component (Talbot & Thrane 2018), dubbed the Power Law + Peak model in Abbott et al. (2021d,c), and that the mass ratio distribution is also a power law, bounded such that the minimum and maximum masses of the secondary component are the same as those of the primary (Fishbach & Holz 2020). The merger rate is allowed to evolve with redshift as a power-law, following the model in Fishbach et al. (2018) and Abbott et al. (2021c). The hyper-parameters governing these mass and redshift distributions are described in Table 1.
Our goal is to determine if there is a correlation between the BBH spin distribution and the distributions of masses or redshifts. To this end, we modify the correlated spin model introduced in Callister et al. (2021) so that the effective aligned spin (Damour 2001;Santamaria et al. 2010), where we use log to indicate log base 10 everywhere and ln to indicate the natural logarithm. The pivot points of z = 0.5 and m 1 = 10 M are chosen to be near the peak of the population distributions of those parameters, but we have verified that the exact values do not affect our results. To avoid unphysically narrow distributions using the model above, we impose a cut such that The full set of Λ χ eff parameters are described in Table 2. We note that unlike the model of Callister et al. (2021), our model in Eq. 2 does not allow for correlations between χ eff and mass ratio. However, we amend the model to incorporate mass ratio correlations later in Section 3. Given a set of posterior samples for the binary parameters θ of N individual BBH events, they can be combined to obtain posteriors on the hyper-parameters governing the population-level distributions: The factor of π(Λ) in Eq. 6 represents the prior probability for the hyper-parameters, given in Tables 1-2. The likelihood in Eq. 7 consists of a Monte Carlo integral over the posterior samples j for each individual event i, where π pop (θ i,j |Λ) is the product of the population distributions for the masses, redshift, and χ eff in Eq. 2. π PE (θ i,j ) is the original prior that was applied during the parameter estimation (PE) of the binary parameters for individual events. Finally, the factor of α(Λ) N accounts for the fact that the observed BBH sources are a biased sample of the underlying astrophysical distribution (Loredo 2004;Mandel et al. 2019;Vitale 2020). This selection effect arises from the dependence of the sensitivity of the detector network on the intrinsic parameters of the source. We use the GWPopulation package ) and the dynesty nested sampler (Speagle 2020) to obtain samples from the hyper-parameter posterior in Eq. 6. We include the 69 BBH events reported in GWTC-3 with a false alarm rate (FAR) less than 1 per year (Abbott et al. 2021b), as was done for the BBH analyses in Abbott et al. (2021c). We use the individual-event posterior samples publicly released by LIGO and Virgo (Abbott et al. 2018(Abbott et al. , 2020a obtained with the IMRPhenomPv2 waveform model for events first published in GWTC-1 (Abbott et al. 2019a), and using a combination of different waveform models including the effects of spin precession and higherorder modes for events reported in later catalogs 1 (Abbott et al. 2021g,a,b). We calculate α(Λ) following the method described in Farr (2019), using the sensitivity estimates for BBH systems released at the end of the most recent LIGO-Virgo observing run (O3) obtained via a simulated injection campaign (Abbott et al. 2021h).

RESULTS FOR GWTC-3
We first allow the χ eff distribution to be correlated with redshift and primary mass individually. We recover mass and redshift hyper-parameter posteriors consistent with those reported in Abbott et al. (2021c). The evolution of µ χ and σ χ as a function of redshift for individual hyper-parameter posterior samples obtained with the model only allowing for redshift correlations is shown in Fig. 1. In this case the δµ m1 and δlog σ m1 parameters are fixed to zero. While µ χ does not exhibit significant evolution with redshift, σ χ increases with redshift. The corresponding posteriors on Λ χ eff are shown in Fig. 2. The δµ z posterior is consistent with 0, indicating no evolution of the mean of the χ eff distribution as a function of redshift, but we find that δlog σ z = 0 is disfavored at 98.6% credibility. We recover δlog σ z = 0.93 +0.54 −0.54 , (maximum posterior value and 90% credible interval calculated with the maximum posterior density method) indicating that the spin distribution broadens with increasing redshift. This correlation can also be visualized by examining slices of the χ eff distribution at fixed redshift, as shown in Fig. 3. The thickness of the 90% credible region for each slice is comparable, but the distribution becomes noticeably broader with increasing redshift. This indicates that the apparent correlation between the width of the spin distribution and the redshift is not dominated by increased uncertainty in the constraint on the distribution at higher redshifts. The location of the peak of the χ eff distribution does not vary significantly between the different redshift slices, consistent with the lack of observed evolution in µ χ .
While these results present compelling evidence for a correlation between χ eff and redshift, it is possible that the BBH mass and redshift distributions are themselves correlated, so that our recovered redshift correlation is actually a manifestation of a mass vs. χ eff correlation viewed through the wrong model. To test this, we now allow the χ eff distribution to be correlated only with primary mass, fixing δµ z = 0, δlog σ z = 0. The posteriors for the spin hyper-parameters for the primary mass correlation model are shown in Fig. 4, and the distributions for χ eff along different slices in primary mass are shown in Fig. 5. We find a similar but less significant trend between the width of the χ eff distribution and the primary mass, where δlog σ m1 = 0.09 +0.10 −0.08 , and δlog σ m1 = 0 is excluded at 93.5% credibility. While the χ eff distributions in Fig. 5 appear visually to increase both in mean and width with increasing primary mass, the posterior on δµ m1 is still consistent with 0 at 31.1% credibility. In Section 4, we comment on the possibility of a correlation between primary mass and χ eff being spuriously introduced due to mismodeling a true correlation between redshift and χ eff . In order to determine if there is a preference in whether the redshift or the primary mass drives the evolution of the width of the χ eff distribution, we now analyze the data with a model that allows for correlations with both parameters. The posteriors on the spin hyper-parameters are shown in Fig. 6. We obtain much weaker constraints on δlog σ z and δlog σ m1 individually, but both posteriors have more support for positive than  Table 3. Natural log Bayes factors comparing the models allowing for correlations between redshift, primary mass, mass ratio, and χ eff and the base model without any correlations. The redshift-only models include both δµz and δlog σz as free parameters and only δlog σz as a free parameter with δµz fixed to zero. The models including a correlation with primary mass are disfavored by the data.
negative values, and the point δlog σ z , δlog σ m1 = (0, 0) is disfavored, lying at the 96.2% credibility contour. Based on these posteriors alone, we cannot confidently identify whether the correlation is dominated by the redshift or the primary mass, but comparing the Bayes factors between the various models that we have considered can provide an indication of which correlation is statistically preferred. Table 3 shows the natural log Bayes factors between the three correlated models we have presented so far and an uncorrelated model, where only µ 0 and log σ 0 in Eq. 2 are left as free parameters. We repeat the analysis allowing only for redshift correlations with fixed δµ z = 0 in order to gauge the effect of the Occam penalty on the Bayes factor, since this sub-hypothesis is consistent with our initial redshift correlation result. While none of the values are particularly statistically significant, we find that the models including a correlation with primary mass are disfavored relative to the redshift-only models. We emphasize that while we have ensured consistent priors between all the models which are sub-hypotheses of each other, the choice of priors for individual hyper-parameters is necessarily somewhat arbitrary, complicating the interpretation of the Bayes factors 2 .
Finally, we want to ensure that the apparent correlation between redshift and χ eff is not falsely introduced by the previously-reported correlation between the mean of the χ eff distribution and mass ratio Abbott et al. 2021c). We modify the model in 2 We can use a simple back-of-the-envelope calculation to determine the effect of the narrower width of the prior on δlog σz on the Bayes factor for the redshift correlation models. Since the prior volume for δlog σz is 2/3.5 times the prior volume for δlog σq and δlog σm 1 , the Bayes factor would be increased for the redshift correlation model by ∼ − ln 2/3.5 ∼ 0.56. This is comparable to the uncertainty in the Bayes factor due to the Monte Carlo integration in the likelihood in Eq. 7.
The priors on δµ q and δlog σ q are given in Table 2, and the Bayes factor between this model and the uncorrlated model is also shown in Table 3. The results we obtain with this model are consistent with both the previouslyreported mass ratio correlation and the redshift correlation presented earlier in this work. The δµ q posterior peaks at negative values, indicating that the mean of the χ eff distribution shifts towards smaller values as the mass ratio becomes more equal. Simultaneously, the δlog σ z posterior peaks at positive values, similar to the posteriors shown in Figs. 2, 6. Thus, we conclude that the increase in the width of the χ eff distribution that we report here is independent of the previously-established correlation between χ eff and mass ratio.

VALIDATION OF RESULTS
The results presented so far suggest that for the 69 BBH events we analyze with FAR < 1 yr −1 detected up to the end of O3, the χ eff distribution broadens with increasing redshift. As noted above, though, it is difficult to disentangle a correlation between spin and redshift from a correlation between spin and mass and to distinguish the extent to which the broadening of the χ eff distribution is due to increased uncertainty in spin measurements at high redshifts. We test the robustness of the observed redshift correlation first with a series of simulated populations, then by analyzing the GWTC-3 data with alternative models allowing for a correlation between χ eff and redshift.

Simulated populations
One potential concern is that such a spin-redshift correlation could be introduced due to the degradation of the constraint on χ eff for individual events at farther redshifts, as these sources are generally detected with smaller signal-to-noise ratios (SNRs). To verify whether such a spurious correlation could be introduced, we simulate a population of 69 BBH events detected by a Hanford-Livingston detector network operating at O3 sensitivity (Abbott et al. 2020b) drawn from an intrinsically uncorrelated population. The true values of the hyper-parameters describing the population are given in Table 4. We consider an event to be "detected" if it has a network optimal SNR ≥ 9. We perform individual-event parameter estimation using the reduced order quadrature implementation (Smith et al. 2016)  . We then recover posteriors on the mass, redshift, and spin hyperparameters using the three population models applied to the real data in the previous section.
The true values of all hyper-parameters describing this simulated population are recovered within at least the 3σ level, although the posteriors for both δlog σ z and δlog σ m1 tend to prefer negative values. This suggests that it is easier to rule out a distribution that broadens with redshift or mass rather than one that becomes narrower. A broadening distribution implies the presence of highly spinning sources at high mass/redshift, which are easier to detect. Therefore, their absence from the observed population indicates that they are also absent from the underlying population. The converse is not true. If the χ eff distribution instead became narrower, there would be an increasing number of sources with nearly zero spin at high mass/redshift. These sources are harder to detect, so their absence from the observed distribution does not necessarily imply their absence from the underlying distribution.
The χ eff distributions along different slices in redshift and primary mass are shown in Fig. 7. These distributions are morphologically distinct from those obtained for the real data in Figs. 3 and 5. The distributions obtained for this uncorrelated population appear to initially get narrower between the lower two values of redshift and primary mass but then broaden again. This indicates that there is no strong evidence for either an increase or a decrease in the width of the χ eff distribution with either parameter. There is more uncertainty in the distributions at higher primary mass and redshift, a feature which is not observed for the real data, further increasing our confidence in our measurement of a broadening in the spin distribution with mass and/or redshift.
The natural log Bayes factors between each of the models allowing for correlations and the model without   Table 5. Natural log Bayes factors comparing the models allowing for correlations between redshift, primary mass, and χ eff and the base model without any correlations for the simulated population with no correlations.
any correlations for this simulated population are given in Table 5. In this case, the redshift-correlation model is disfavored at a level comparable to the primary-masscorrelation model, as expected since the true population does not contain any correlations. This simulation suggests that a positive correlation between the width of the χ eff distribution and either the redshift or primary mass is unlikely to be spuriously introduced by the worsening constraint on χ eff for individual events at higher redshift or primary mass. We next seek to verify if we are able to successfully recover a correlation in a simulated population similar to the one we find in real data, and whether a redshift correlation can manifest as a primary mass correlation when analyzed with the wrong model. We again generate 69 sources detected at O3 sensitivity now drawn from a distribution with δlog σ z = 0.85. All the other true hyper-parameter values are the same as for the uncorrelated population given in Table 4. We initially analyze this simulated population allowing for mass and redshift correlations in turn. As before, all the hyper-parameter values are recovered within 3σ credibility, although the posterior on δlog σ z peaks below the true value. The distributions for χ eff along different slices in mass and redshift shown in Fig. 8 are similar to those recovered in the real data. They become consistently wider between increasing values of redshift and primary mass, although they also become more uncertain, which is not observed in the real data. This is because the posteriors for individual events in the simulated data set only extend up to z = 1.14, m 1 = 104 M , while those for real events include values up to z = 1.90, m 1 = 296 M , meaning there is more resolving power for the mass and redshift distributions at large masses and redshifts in the real data. While the posterior on δlog σ m1 includes 0 within the 81.4% credible level, this simulation indicates that a true correlation between redshift and χ eff could be perceived as a correlation between primary mass and χ eff if analyzed with the wrong model. It also reinforces the significance of our finding a preference for positive values of δlog σ z in the real data, since both simulations tend to recover smaller values of δlog σ z compared to the true value. This is further supported by the Bayes factors we recover for the correlated simulation, given in Table 6. In this case, the recovered correlation with redshift-only is not significant enough to overcome the Occam penalty for adding parameters relative to the uncorrelated model. However, the relative Bayes factors between the three models considered are similar to those obtained for the real data. This indicates that the redshift-only correlation found in the real data is preferred over the other correlated models we explore at a similar level of statistical significance to the simulated population with a known correlation in δlog σ z . On the other hand, the relative Bayes factors for the uncorrelated simulation span a much narrower range, revealing  Table 6. Natural log Bayes factors comparing the models allowing for correlations between redshift, primary mass, and χ eff and the base model without any correlations for the simulated population with δlog σz = 0.85.
that none of the correlated models is strongly preferred over the others in the case of no correlation. Conversely, it is also possible that a true correlation between primary mass and χ eff could manifest as a correlation between redshift and χ eff if analyzed with the wrong model. To verify whether this is a potential false source of the redshift-χ eff correlation we observe, we simulate a third population of 69 events detected at O3 sensitivity with hyper-parameters identical to the previous two, except now δlog σ z = 0 and δlog σ m1 = 0.15. We find that when this population is analyzed with the model allowing for a redshift correlation only, δlog σ z = 0 is incorrectly ruled out at 99.0% credibility, indicating that a primary mass-χ eff correlation can indeed be confused for a redshift-χ eff correlation. This is similar to the case explored in the previous two paragraphs where δlog σ m1 = 0 was disfavored at the 81.4% credible level for a true correlation between redshift and spin analyzed with the wrong model. However, we find in general that a correlation between spin and primary mass is easier to confidently identify with 69 events detected at O3 sensitivity. The posterior on δlog σ m1 disfavors 0 at 99.7% credibility under the model allowing only for a primary mass correlation and at 99.5% credibility under the model allowing for both primary mass and redshift correlations. This is illustrated in the corner plot of the posteriors on δlog σ m1 , δlog σ z in Fig. 9 for both the mass-correlated population and the redshift-correlated population previously presented in Fig. 8. While it is difficult to experimentally distinguish between mass and redshift evolution of the spin distribution in the case of a redshiftcorrelated population, the same ambiguity is not present for the mass-correlated population. This further supports the redshift-χ eff correlation we find in the real data, since a correlation between primary mass and χ eff would likely have been unambiguously identified.

Alternative models
While the results of the two simulated populations lend credence to the finding that the width of the χ eff distribution increases with redshift in the real data, we want to verify whether this result is model-driven. To this end, we now analyze the GWTC-3 events with a different model-a mixture of truncated Gaussians with a redshift-dependent mixing fraction, where We choose three different priors on the µ and σ parameters for each Gaussian, outlined in Table 7, as proxies for different astrophysical scenarios. When modeling the effective spin distribution with a single, redshift-independent Gaussian, it is found that the mean effective spin is χ eff = 0.06 (Miller et al. 2020;Abbott et al. 2021c). For our first prior choice, we assume that this result encapsulates the bulk of the population, fixing µ b = 0.06, and let the second mixture component model the departure from this assumption  Table 7. The solid black line shows the mean, while the dashed black lines show the 90% credible region.
as we look to higher redshifts. We allow this secondary sub-population to peak at either positive or negative values of χ eff , but the absolute value of the peak must be ≥ 0.1, so as to minimize the degeneracy between the bulk and this second sub-population. The evolution of the mixture fraction with redshift for individual hyper-parameter posterior samples is shown in Fig. 10. The posterior on the mixture fraction differs substantially from the prior, which is symmetric about f = 0.5 for all redshifts. This result indicates that the contribution of the Gaussian defined by the parameters µ a and σ a becomes more significant as redshift increases.
The posterior for µ a is 3.4 times more likely to be positive than negative, and the positive part of the posterior rails against the boundary at µ a = 0.1, indicating that it may be trying to replicate the bulk distribution rather than to extract an independent component. σ a is only weakly constrained to be > −1.02 at 90% credibility. Nonetheless, the χ eff distributions at different redshift slices shown in the top of Fig. 11 for this model and prior choice paint a similar picture to the linear evolution model presented in Section 3; although skewed towards positive values of χ eff , the distributions become wider with increasing redshift. The second prior choice targets two populations with χ eff distributions characterized by different widths. The peak of the narrower bulk distribution is now fixed to µ b = 0 to capture the BBH population formed dynamically with small spin and the majority of field binaries predicted to have negligible spin (Fuller & Ma 2019). Meanwhile the posterior on the peak of the broader Gaussian, µ a , is largely uninformative. The posteriors for the width parameters and the mixture fraction are very similar to those obtained with the first prior, with  σ a > −0.25 at 90% credibility and σ b more strongly constrained to peak at −1.23 +0.17 −0.17 . These results do not provide significant evidence for two distinct populations, but the χ eff distributions at different redshift slices shown in the middle panel of Fig. 11 still exhibit some broadening with increasing redshift.
The final prior choice is designed to look for a redshiftdependent excess of preferentially-aligned BBH systems. The mean of the bulk distribution is still fixed to µ b = 0, as expected for a population of dynamically-formed systems with isotropic spin tilts (Portegies Zwart & McMillan 2002;Antonini & Rasio 2016;Rodriguez et al. 2015;Gerosa & Berti 2017;Rodriguez et al. 2019;Kimball et al. 2020;Gerosa & Fishbach 2021). The mean of the aligned population is restricted to 0.1 < µ a < 1, in order to represent systems formed via isolated binary evolution in the field (Gerosa et al. 2018;Qin et al. 2018;Zaldarriaga et al. 2018;Belczynski et al. 2020;Bavera et al. 2020Bavera et al. , 2021Bavera et al. , 2022, and the corresponding Gaussian is truncated on [0, 1] instead of [−1, 1]. The χ eff distributions at different redshift slices for this prior choice are shown in the bottom panel of Fig. 11. The population model is asymmetric by construction, so the mean of the distribution also increases with redshift along with the width.
The posterior on the mixture fraction is again very similar to the one obtained with the first prior choice in Fig. 10. Since the mixture fraction increases with redshift for all three prior choices, the second Gaussian added on top of the bulk distribution always contributes more at high redshifts, consistent with the conclusion that the χ eff distribution broadens with increasing redshift. These results indicate that the apparent correlation between the width of the χ eff distribution and the redshift is not driven by our choice of linear evolution model in Section 2, but can be observed under a variety of population models and priors.

DISCUSSION AND CONCLUSION
In this work, we have found weak but robust evidence for an increase in the width of the χ eff distribution of BBH systems with increasing redshift. We have verified that this correlation is unlikely to be spuriously introduced by the increased uncertainty in the χ eff posteriors for individual high-redshift sources using a simulated population and that this trend remains present using alternative population models. We also observe a less significant correlation between the width of the χ eff distribution and the primary mass, although we find that such a correlation can be falsely recovered if a population with a χ eff -redshift correlation is instead analyzed assuming only a χ eff -primary mass correlation. When allowing for correlations with both redshift and primary mass, we find evidence for an increase in the width of the χ eff distribution with one or the other but cannot distinguish which.
A correlation between χ eff and redshift might be explained by one of two broad scenarios. First, the correlation could be indicative of multiple sub-populations arising from distinct formation channels, each of which occupy a different region in the (χ eff , z) plane. The third model defined in Section 4.2, for example, serves as a proxy for this scenario. In this case, the symmetric χ eff component centered at zero might serve to capture binaries assembled dynamically in dense stellar environments (Portegies Zwart & McMillan 2002;Antonini & Rasio 2016;Rodriguez et al. 2015;Gerosa & Berti 2017;Rodriguez et al. 2019;Kimball et al. 2020;Gerosa & Fishbach 2021), while the preferentially-positive χ eff component corresponds to a sub-population of systems formed via isolated binary evolution in the field (Gerosa et al. 2018;Qin et al. 2018;Zaldarriaga et al. 2018;Belczynski et al. 2020;Bavera et al. 2020Bavera et al. , 2021Bavera et al. , 2022. However, such a mixture between two sub-populations generally leads to a shift in the mean of the χ eff distribution with redshift. Under the more generic model explored in Section 3, we much more strongly prefer evolution of the width with redshift, although evolution of µ χ is not strictly ruled out. The second possibility, broadly, is that a spin-redshift correlation arises due to evolutionary processes operating within a single population. Within isolated binary evolution, for example, tidal interactions may be responsible for correlating black hole spins and redshifts. Some  Table 7 (1-3 from top to bottom).
authors predict that isolated black holes are naturally born slowly rotating due to efficient angular momentum transport from stellar cores (Spruit 2002;Fuller & Ma 2019). Spin can nevertheless be introduced via tidal torques exerted by the first-born black hole on its companion, prior to the companion's core collapse (Zaldarriaga et al. 2018;Qin et al. 2018;Bavera et al. 2021Bavera et al. , 2020Fuller & Lu 2022). In this scenario, the shortest-period binaries both acquire the largest spins and merge most promptly, correlating the spins of binary black holes and the redshifts at which they merge (Qin et al. 2018;Fuller & Ma 2019). This effect may be enhanced by the lower metallicities present at high redshifts, which diminish the strength of stellar winds and hence prevent binary orbits from widening and avoiding tidal spin-up. Bavera et al. (2022) corroborate this prediction using population synthesis simulations, finding that an increasing fraction of systems undergo tidal spin-up in close orbits at high redshift.
Because there are systems that do not meet the criteria for efficient tidal spin-up at all redshifts, they find that the spin distribution both broadens and increases in mean. Another possible explanation comes from the effect of metallicity on the efficiency of angular momentum transport in the envelope of the stellar progenitors of the BBH system. Because stellar winds are weaker at low metallicity (higher redshift), this leads to a less efficient removal of the angular momentum stored in the stellar envelope and hence a possibly higher spin of the resulting black hole (Qin et al. 2018;Fuller & Ma 2019). However, this naive picture is complicated by the interplay between the strength of stellar winds, the extent of the hydrogen-burning region inside the star, and the efficiency of elemental mixing within the star. As such, Belczynski et al. (2020) find a non-monotonic relationship between the metallicity and the black hole spin for a given carbon-oxygen-core mass of the progenitor. This relationship further depends on the stellar evolution model employed in the binary evolution simulation. Because of the uncertainty in the processes that impart natal spin to the black hole, we cannot place meaningful constraints on these models.
While both of the possibilities described above lead to systems with larger spin magnitudes at increasing redshift, we must also explain the increased number of systems with negative values of χ eff due to the symmetric broadening of the distribution. One potential explanation within the framework of isolated binary evolution is that natal kicks are stronger at higher redshift, leading to more significant misalignment of the BH spin relative to the orbital angular momentum upon birth. It is not clear how such an effect would arise, however. Furthermore, even with moderately strong (∼ 100km s −1 ) natal kicks, the median of the χ eff distribution for systems formed in the field is still much higher than what we recover (Gerosa et al. 2018). Alternatively, the mixing fraction between systems formed via isolated evolution and dynamical assembly may remain constant with redshift, but the spin magnitudes of systems formed via both channels can increase due to one or both of the scenarios previously discussed. In this case, the orientations of the spins in systems formed dynamically would remain isotropic, but their magnitudes would be larger at higher redshift, leading to a broadening of the distribution but not a shift in the mean.
One caveat with our analysis is that we restrict any correlations to be with χ eff , rather than also allowing for mass-redshift correlations. It is likely that the complex processes leading to the formation of BBH systems would produce correlations between all of these parameters, rather than just with spin (Qin et al. 2018;Bavera et al. 2021;Fuller & Lu 2022;. However, it is difficult to encapsulate all possible correlations in a simple phenomenological model without dramatically increasing the dimensionality of the parameter space. We leave the exploration of simultaneous correlations between all parameters to future work.