The Binary Black Hole Spin Distribution Likely Broadens with Redshift

The growing catalog of compact binary mergers detected with gravitational waves (Abbott et al. 2019a, 2021a, 2021b) 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; Bouffanais et al. 2021; Zevin et al. 2021). The latest data from the Advanced LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) observatories have revealed that the BBH mass distribution has substructure (Li et al. 2021; Tiwari & Fairhurst 2021; Veske et al. 2021; Edelman et al. 2022) beyond a smooth power law and single peak or break at ∼40 M_⊙ (Fishbach & Holz 2017; Talbot & Thrane 2018; Abbott et al. 2019b). BBH spins are found to be small but nonzero (Roulet & Zaldarriaga 2019; Wysocki et al. 2019; Miller et al. 2020; Biscoveanu et al. 2021; García-Bellido et al. 2021), with some of their tilts misaligned to the orbital angular momentum (Talbot & Thrane 2017; Abbott et al. 2019b, 2021c), although these conclusions have been challenged when applying a different population model (Galaudage et al. 2021; Roulet 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; Abbott et al. 2021c; Callister et al. 2021; Franciolini & Pani 2022; Tiwari 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; Miller & Lauburg 2009; McKernan et al. 2012; Benacquista & Downing 2013; Rodriguez et al. 2015; Antonini & Rasio 2016; Fishbach et al. 2017; 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 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 (Bavera et al. 2021; Zevin & Bavera 2022), so may hint at the superposition of multiple populations formed via independent channels with unique q versus χ_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 (${ \mathcal M }={\left({m}_{1}{m}_{2}\right)}^{3/5}/{\left({m}_{1}+{m}_{2}\right)}^{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, 2020b, 2020c; Franciolini et al. 2022), 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; Fishbach et al. 2021). Heavier black holes are predicted to form from lower-metallicity stellar progenitors at higher redshifts, which could lead to such a correlation (Mapelli et al. 2019; Neijssel et al. 2019; Belczynski et al. 2020). 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. 2022). Motivated by these theoretical predictions, Fishbach et al. (2021) find that redshift dependence 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 that 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.

We employ the framework of hierarchical Bayesian inference to constrain the hyperparameters 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₁, 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. (2021c, 2021d), 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 hyperparameters 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; Ajith et al. 2011; Ajith 2011),

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{{\chi }_{1}\cos {\theta }_{1}+q{\chi }_{2}\cos {\theta }_{2}}{1+q},\end{eqnarray} \tag{ 1 }$$

is modeled as a truncated Gaussian on [−1, 1] whose mean and variance can each evolve linearly with primary mass and redshift,

$$\begin{eqnarray}&&{\pi }_{\mathrm{pop}}({\chi }_{\mathrm{eff}}| {{\rm{\Lambda }}}_{{\chi }_{\mathrm{eff}}},{m}_{1},z)={ \mathcal N }({\chi }_{\mathrm{eff}};{\mu }_{\chi },{\sigma }_{\chi }),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\mu }_{\chi }(z,{m}_{1})={\mu }_{0}+\delta {\mu }_{z}(z-0.5)+\delta {\mu }_{{m}_{1}}\left(\displaystyle \frac{{m}_{1}}{10\,{M}_{\odot }}-1\right)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{log}}{{\sigma }}_{{\chi }}(z,{m}_{1}) & = & {\rm{log}}{{\sigma }}_{0}+{\delta }{\rm{log}}{{\sigma }}_{z}(z-0.5)\\ & & +{\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}\left(\displaystyle \frac{{m}_{1}}{10\,{M}_{\odot }}-1\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where we use $\mathrm{log}$ to indicate log base 10 everywhere and $\mathrm{ln}$ to indicate the natural logarithm. The pivot points of z = 0.5 and m₁ = 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

$$\begin{eqnarray}&&{\pi }_{\mathrm{pop}}({\chi }_{\mathrm{eff}}| {{\rm{\Lambda }}}_{{\chi }_{\mathrm{eff}}})=\left\{\begin{array}{l}0,\quad \mathrm{if}\,{\mathrm{log}}_{10}{\sigma }_{\chi }(z,{m}_{1})\lt -3\\ { \mathcal N }({\chi }_{\mathrm{eff}};{\mu }_{\chi },{\sigma }_{\chi }),\quad \mathrm{elsewhere}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

The full set of ${{\rm{\Lambda }}}_{{\chi }_{\mathrm{eff}}}$ parameters are described in Table 2. We note that unlike the model of Callister et al. (2021), our model in Equation (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.

Table 2. Spin Hyperparameter Priors and Descriptions for the Model in Equation (2)

Parameter	Description	Prior
μ0	Independent χeff mean	U(−1, 1)
$\mathrm{log}{\sigma }_{0}$	Log of independent χeff width	U(−1.5, 0.5)
δ μq	q-dependent χeff mean	U(−2.5, 1)
${\delta }{\rm{log}}{{\sigma }}_{q}$	Log of q-dependent χeff width	U(−2, 1.5)
δ μz	z-dependent χeff mean	U(−2.5, 1)
${\delta }{\rm{log}}{{\sigma }}_{z}$	Log of z-dependent χeff width	U(−0.5, 1.5)
$\delta {\mu }_{{m}_{1}}$	m1-dependent χeff mean	U(−2.5, 1)
${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$	Log of m1-dependent χeff width	U(−2, 1.5)

Download table as:  ASCIITypeset image

Given a set of posterior samples for the binary parameters θ of N individual BBH events, they can be combined to obtain posteriors on the hyperparameters governing the population-level distributions:

$$\begin{eqnarray}&&p({\rm{\Lambda }}| \{d\})\propto { \mathcal L }(\{d\}| {\rm{\Lambda }})\pi ({\rm{\Lambda }}),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{ \mathcal L }(\{d\}| {\rm{\Lambda }})\propto \displaystyle \frac{1}{{\alpha }{\left({\rm{\Lambda }}\right)}^{N}}\prod _{i}^{N}\sum _{j}\displaystyle \frac{{{\pi }}_{{\rm{p}}{\rm{o}}{\rm{p}}}({{\theta }}_{i,j}| {\rm{\Lambda }})}{{{\pi }}_{{\rm{P}}{\rm{E}}}({{\theta }}_{i,j})}.\end{eqnarray} \tag{ 7 }$$

The factor of π(Λ) in Equation (6) represents the prior probability for the hyperparameters, given in Tables 1–2. The likelihood in Equation (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 Equation (2); π_PE(θ_i,j ) is the original prior that was applied during the parameter estimation of the binary parameters for individual events. Finally, the factor of $\alpha {\left({\rm{\Lambda }}\right)}^{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; Thrane & Talbot 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 (Talbot et al. 2019) and the dynesty nested sampler (Speagle 2020) to obtain samples from the hyperparameter posterior in Equation (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, 2020a, 2021e, 2021f) 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 higher-order modes for events reported in later catalogs ⁵ (Abbott et al. 2021a, 2021b, 2021g). 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).

We first allow the χ_eff distribution to be correlated with redshift and primary mass individually. We recover mass and redshift hyperparameter posteriors consistent with those reported in Abbott et al. (2021c). The evolution of μ_χ and σ_χ as a function of redshift for individual hyperparameter posterior samples obtained with the model only allowing for redshift correlations is shown in Figure 1. In this case the $\delta {\mu }_{{m}_{1}}$ and ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$ parameters are fixed to zero. While μ_χ does not exhibit significant evolution with redshift, σ_χ increases with redshift. The corresponding posteriors on ${{\rm{\Lambda }}}_{{\chi }_{\mathrm{eff}}}$ are shown in Figure 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 ${\delta }{\rm{log}}{{\sigma }}_{z}=0$ is disfavored at 98.6% credibility. We recover ${\delta }{\rm{log}}{{\sigma }}_{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 Figure 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 μ_χ .

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 versus χ_eff correlation viewed through the wrong model. To test this, we now allow the χ_eff distribution to be correlated only with primary mass, fixing ${\delta }{{\mu }}_{z}=0,{\delta }{\rm{log}}{{\sigma }}_{z}=0$. The posteriors for the spin hyperparameters for the primary mass correlation model are shown in Figure 4, and the distributions for χ_eff along different slices in primary mass are shown in Figure 5. We find a similar but less significant trend between the width of the χ_eff distribution and the primary mass, where ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}={0.09}_{-0.08}^{+0.10}$, and ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}=0$ is excluded at 93.5% credibility. While the χ_eff distributions in Figure 5 appear visually to increase both in mean and width with increasing primary mass, the posterior on $\delta {\mu }_{{m}_{1}}$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 hyperparameters are shown in Figure 6. We obtain much weaker constraints on ${\delta }{\rm{log}}{{\sigma }}_{z}$ and ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$ individually, but both posteriors have more support for positive than negative values, and the point ${\delta }{\rm{log}}{{\sigma }}_{z},{\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}=(0,0)$ is disfavored, lying at the 96.2% credibility contour.

Zoom In Zoom Out Reset image size

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 μ₀ and $\mathrm{log}{\sigma }_{0}$ in Equation (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 subhypothesis 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 subhypotheses of each other, the choice of priors for individual hyperparameters is necessarily somewhat arbitrary, complicating the interpretation of the Bayes factors. ⁶

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

Correlation	$\mathrm{ln}(\mathrm{BF})$
Redshift, δ μz = 0	1.43
Redshift	0.24
Primary mass	−2.09
Mass ratio	3.75
Redshift & primary mass	−3.63
Redshift & mass ratio	3.20

Note. The redshift-only models include both δ μ_z and ${\delta }{\rm{log}}{{\sigma }}_{z}$ as free parameters and only ${\delta }{\rm{log}}{{\sigma }}_{z}$ as a free parameter with δ μ_z fixed to zero. The models including a correlation with primary mass are disfavored by the data.

Download table as:  ASCIITypeset image

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; Callister et al. 2021). We modify the model in Equation (2) to allow for correlations with both redshift and mass ratio rather than redshift and primary mass:

$$\begin{eqnarray}&&{\mu }_{\chi }(z,q)={\mu }_{0}+\delta {\mu }_{z}(z-0.5)+\delta {\mu }_{q}\left(q-0.7\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{log}}}_{10}{{\sigma }}_{{\chi }}(z,q) & = & {\rm{log}}{{\sigma }}_{0}+{\delta }{\rm{log}}{{\sigma }}_{z}(z-0.5)\\ & & +{\delta }{\rm{log}}{{\sigma }}_{q}\left(q-0.7\right).\end{array}\end{eqnarray} \tag{ 9 }$$

The priors on δ μ_q and ${\delta }{\rm{log}}{{\sigma }}_{q}$ are given in Table 2, and the Bayes factor between this model and the uncorrelated model is also shown in Table 3. The results we obtain with this model are consistent with both the previously reported 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 toward smaller values as the mass ratio becomes more equal. Simultaneously, the ${\delta }{\rm{log}}{{\sigma }}_{z}$ posterior peaks at positive values, similar to the posteriors shown in Figures 2 and 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.

The results presented so far suggest that for the 69 BBH events we analyze with FAR < 1 yr⁻¹ 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.

4.1. 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 hyperparameters 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) of the IMRPhenomPv2 waveform (Hannam et al. 2014; Husa et al. 2016; Khan et al. 2016) and the dynesty nested sampler (Speagle 2020) through the bilby package (Ashton et al. 2019; Romero-Shaw et al. 2020). 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.

Table 4. True and Recovered Values of the Hyperparameters Describing the Mass, Spin, and Redshift Distributions from which the Simulated Population with No Correlation was Drawn

Parameter	Value	Recovery
α	4	${4.74}_{-0.98}^{+1.33}$
β	1.5	${1.02}_{-1.44}^{+2.10}$
mmax	50 M⊙	${46.63}_{-7.87}^{+46.35}\,{M}_{\odot }$
mmin	5 M⊙	${4.99}_{-1.10}^{+0.87}\,{M}_{\odot }$
δm	5 M⊙	${4.88}_{-2.12}^{+4.12}\,{M}_{\odot }$
μm	35 M⊙	${31.40}_{-1.59}^{+1.59}\,{M}_{\odot }$
σm	4 M⊙	${3.09}_{-2.09}^{+0.99}\,{M}_{\odot }$
λ	0.04	${0.012}_{-0.012}^{+0.019}$
λz	3	${5.15}_{-1.83}^{+1.83}$
μ0	0.05	${0.025}_{-0.088}^{+0.076}$
$\mathrm{log}{\sigma }_{0}$	−0.85	$-{0.82}_{-0.20}^{+0.22}$
δ μz	0	$-{0.17}_{-0.11}^{+0.21}$
${\delta }{\rm{log}}{{\sigma }}_{z}$	0	$-{0.24}_{-0.47}^{+0.62}$
$\delta {\mu }_{{m}_{1}}$	0	${0.023}_{-0.035}^{+0.035}$
${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$	0	$-{0.14}_{-0.08}^{+0.07}$

Note. The correlated population is described by the same hyperparameters with the exception of ${\delta }{\rm{log}}{{\sigma }}_{z}=0.85$. The recovered values are represented by the maximum posterior value and 90% credible interval calculated with the maximum posterior density method obtained with the model that allows for both primary mass and redshift correlations

Download table as:  ASCIITypeset image

The true values of all hyperparameters describing this simulated population are recovered within at least the 3σ level, although the posteriors for both ${\delta }{\rm{log}}{{\sigma }}_{z}$ and ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$ 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 Figure 7. These distributions are morphologically distinct from those obtained for the real data in Figures 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 that 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.

Zoom In Zoom Out Reset image size

The natural log Bayes factors between each of the models allowing for correlations and the model without 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-mass-correlation 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.

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

Correlation	$\mathrm{ln}(\mathrm{BF})$
Redshift	−1.93
Primary mass	−1.59
Both	−3.35

Download table as:  ASCIITypeset image

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 ${\delta }{\rm{log}}{{\sigma }}_{z}=0.85$. All the other true hyperparameter 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 hyperparameter values are recovered within 3σ credibility, although the posterior on ${\delta }{\rm{log}}{{\sigma }}_{z}$ peaks below the true value. The distributions for χ_eff along different slices in mass and redshift shown in Figure 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₁ = 104 M_⊙, while those for real events include values up to z = 1.90, m₁ = 296M_⊙, meaning there is more resolving power for the mass and redshift distributions at large masses and redshifts in the real data.

Zoom In Zoom Out Reset image size

While the posterior on ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$ 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 ${\delta }{\rm{log}}{{\sigma }}_{z}$ in the real data, since both simulations tend to recover smaller values of ${\delta }{\rm{log}}{{\sigma }}_{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 ${\delta }{\rm{log}}{{\sigma }}_{z}$. On the other hand, the relative Bayes factors for the uncorrelated simulation span a much narrower range, revealing that none of the correlated models are strongly preferred over the others in the case of no correlation.

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 ${\delta }{\rm{log}}{{\sigma }}_{z}=0.85$

Correlation	$\mathrm{ln}(\mathrm{BF})$
Redshift	−2.45
Primary mass	−4.55
Both	−6.97

Download table as:  ASCIITypeset image

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 hyperparameters identical to the previous two, except now ${\delta }{\rm{log}}{{\sigma }}_{z}=0$ and ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}=0.15$. We find that when this population is analyzed with the model allowing for a redshift correlation only, ${\delta }{\rm{log}}{{\sigma }}_{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 ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}=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 ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}}$ 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 ${\delta }{\rm{log}}{{\sigma }}_{{m}_{1}},{\delta }{\rm{log}}{{\sigma }}_{z}$ in Figure 9 for both the mass-correlated population and the redshift-correlated population previously presented in Figure 8. While it is difficult to experimentally distinguish between mass and redshift evolution of the spin distribution in the case of a redshift-correlated 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.

Zoom In Zoom Out Reset image size

4.2. 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,

$$\begin{eqnarray}&&\pi ({\chi }_{\mathrm{eff}}| {{\rm{\Lambda }}}_{{\chi }_{\mathrm{eff}}})={ \mathcal N }({\chi }_{\mathrm{eff}};{\mu }_{b},{10}^{{\sigma }_{b}})(1-f(z))\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&+\,{ \mathcal N }({\chi }_{\mathrm{eff}};{\mu }_{a},{10}^{{\sigma }_{a}})f(z),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&f(z)=\displaystyle \frac{1}{1+\exp \left[{f}_{A}+{f}_{B}(z-0.5)\right]}.\end{eqnarray} \tag{ 12 }$$

We choose three different priors on the μ and σ parameters for each Gaussian, outlined in Table 7, as proxies for different astrophysical scenarios.

Table 7. Spin Hyperparameter Priors and Descriptions for the Model in Equation (12)

Parameter	Description	Prior 1	Prior 2	Prior 3
μa	Mean of the secondary Gaussian	U(−1, −0.1) ∪ U(0.1, 1)	U(−1, 1)	U(0.1, 1)
μb	Mean of the bulk Gaussian	0.06	0	0
σa	Log-width of secondary Gaussian	U(−1.5, 0.5)	U(−0.7, 1.5)	U(−1.5, 0.5)
σb	Log-width of the bulk Gaussian	U(−1.5, 0.5)	U(−2, −0.7)	U(−1.5, 0.5)
${\chi }_{a,\min }$	Lower bound of the secondary Gaussian	−1	−1	0
fA	Independent offset in mixing fraction	${ \mathcal N }(0,1.5)$	${ \mathcal N }(0,1.5)$	${ \mathcal N }(0,1.5)$
fB	z-dependent offset in mixing fraction	${ \mathcal N }(0,1.5)$	${ \mathcal N }(0,1.5)$	${ \mathcal N }(0,1.5)$

Download table as:  ASCIITypeset image

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 as we look to higher redshifts. We allow this secondary subpopulation 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 subpopulation. The evolution of the mixture fraction with redshift for individual hyperparameter posterior samples is shown in Figure 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 Figure 11 for this model and prior choice paint a similar picture to the linear evolution model presented in Section 3; although skewed toward positive values of χ_eff, the distributions become wider with increasing redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 Figure 11 still exhibit some broadening with increasing redshift.

The final prior choice is designed to look for a redshift-dependent 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; Rodriguez et al. 2015; Antonini & Rasio 2016; 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; Bavera et al. 2020; Belczynski et al. 2020; Bavera et al. 2021, 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 Figure 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 Figure 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.

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 an χ_eff-redshift correlation is instead analyzed assuming only an χ_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 subpopulations 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; Rodriguez et al. 2015; Antonini & Rasio 2016; Gerosa & Berti 2017; Rodriguez et al. 2019; Kimball et al. 2020; Gerosa & Fishbach 2021), while the preferentially positive χ_eff component corresponds to a subpopulation of systems formed via isolated binary evolution in the field (Gerosa et al. 2018; Qin et al. 2018; Zaldarriaga et al. 2018; Bavera et al. 2020; Belczynski et al. 2020; Bavera et al. 2021, 2022). However, such a mixture between two subpopulations 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 authors predict that isolated black holes are naturally born slowly rotating due to efficient angular momentum transport from stellar cores (Spruit 2002; Fuller et al. 2019; 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 (Qin et al. 2018; Zaldarriaga et al. 2018; Bavera et al. 2020, 2021; Fuller & 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 redshifts. 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 nonmonotonic 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 (∼100 km s⁻¹) 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 redshifts, 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; Zevin & Bavera 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.

The authors thank Thomas Dent and Colm Talbot for useful discussions and Michael Zevin, Maya Fishbach, and José Nuño for helpful comments on the manuscript. S.B., C-J.H., K.K.Y.N., and S.V. acknowledge support of the National Science Foundation and the LIGO Laboratory.

LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0757058. This research has made use of data, software, and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. S.V. is supported by the NSF through award PHY-2045740. S.B. is also supported by the NSF Graduate Research Fellowship under grant No. DGE-1122374. The Flatiron Institute is a division of the Simons Foundation, supported through the generosity of Marilyn and Jim Simons. The authors are grateful for computational resources provided by the LIGO Lab and supported by NSF grants PHY-0757058 and PHY-0823459. This paper carries LIGO document number LIGO-P2200105.