Gravitational-wave observations of binary black holes have opened a new window onto massive-star evolution1,2,3,4,5, but population inferences remain hampered by uncertainties in binary physics and initial conditions (see, for example, refs. 6,7,8,9). A central issue is whether (pulsational) pair-instability supernovae (PISNs) carve out a gap in the black hole birth mass distribution (see, for example, refs. 10,11,12); theory predicts pulsations for He cores of ~40–65 M⊙ and full disruption above ~65 M⊙, suppressing black hole formation in the ~40–130-M⊙ range13,14,15,16,17,18. Gravitational-wave observations have so far revealed no sharp deficit of black holes in this mass range (the so-called PISN mass gap)2,3,19,20,21,22,23, motivating scenarios that populate the gap (see, for example, refs. 11,15,24,25), or raising the possibility that the gap may not exist at all (see, for example, ref. 26).
Dynamical environments (for example, globular and nuclear clusters or active galactic nucleus disks) can produce merger remnants that merge again, yielding higher spins27,28,29,30 with isotropic orientations31,32,33, and populating the PISN mass gap. For binaries in which the primary component was produced by a previous merger, the effective combination of the two component spins projected parallel to the orbital angular momentum34,35, χeff (the best measured spin parameter from data), is expected to be broad and symmetric around zero. Reference 36 showed that this distribution is largely independent of model assumptions and uncertainties, and derived an approximately uniform form with ∣χeff∣ ≲ 0.5. This bound follows from the fact that merger remnants are expected to have a nearly universal dimensionless spin magnitude, arem ≃ 0.7, as predicted by general relativity27. Assuming a negligibly spinning companion, χeff satisfies \(| {\chi }_{{\rm{eff}}}| \lesssim \frac{{m}_{{\rm{rem}}}{a}_{{\rm{rem}}}}{{m}_{{\rm{rem}}}+{m}_{2}},\) where mrem and m2 are the remnant and secondary black hole masses, respectively. For m2 ≃ 0.5 mrem, as expected in dynamical formation37, this yields ∣χeff∣ ≲ 0.47. The presence and location of the PISN lower mass limit \(\widetilde{m}\) can therefore be inferred from the primary mass m1 at which the χeff distribution transitions to this broad, symmetric form.
We perform hierarchical Gaussian-process population inference on the fourth LIGO–Virgo–KAGRA gravitational-wave transient catalogue (GWTC-4)38,39 to map black hole spin as a function of primary mass. With the source catalogue now more than twice as large, we obtain tight constraints and are able to probe new features of the population. We fit the χeff distribution to a mixture model comprising a Gaussian distribution, representing the bulk of the population at \({m}_{1}\lesssim \widetilde{m}\), and a higher-mass spin distribution described via a non-parametric Gaussian-process prior.
We identify a transition at \(\widetilde{m}=45.{6}_{-5.6}^{+12.7}\,{M}_{\odot }\) (90% confidence), separating the two populations. The differential merger rates as a function of primary mass for the two populations are given in Fig. 1, while the inferred spin distributions are given in Supplementary Information. Below \(\widetilde{m}\), the data are well described by a single narrow Gaussian χeff distribution, \({\log }_{10}\,\sigma =-1.1{5}_{-0.15}^{+0.13}\), with small and positive mean, \(\mu =0.0{4}_{-0.02}^{+0.02}\), consistent with first-generation black holes. The merger rate of this population drops to zero above \(\widetilde{m}\), implying a gap in the mass spectrum of first-generation mergers. In contrast, the population above \(\widetilde{m}\) exhibits a much broader spin distribution with a median value \(\langle {\chi }_{{\rm{eff}}}\rangle =0.1{1}_{-0.243}^{+0.197}\). Under the proposed model, the data are consistent with a χeff distribution that is symmetric about zero and with the expected distribution of second-generation mergers formed dynamically in dense stellar environments.
Fig. 1: The primary black hole mass spectrum.
The alternative text for this image may have been generated using AI.
Merger rate as a function of primary black hole mass in binaries below (red) and above (grey) the truncation mass \(\widetilde{m}\) separating low- and high-spin populations, calculated at a reference redshift z = 0.2. Given the merger-rate density per unit primary mass at the reference redshift, Γ(m1), the primary-mass distribution for the two populations below and above \(\widetilde{m}\) are reconstructed from the posterior samples as \({\varGamma }_{ < \widetilde{m}}({m}_{1})=\left[1-\eta ({m}_{1})\right]\varGamma ({m}_{1})\) and \({\varGamma }_{ > \widetilde{m}}({m}_{1})=\eta ({m}_{1})\varGamma ({m}_{1})\), respectively. Here η(m1) is a sigmoid mixing function that sets the relative contribution of the two population components as a function of primary mass, with \(\eta (\widetilde{m})=0.9\). This model yields a total merger rate at this redshift of \(33.{4}_{-8.4}^{+13.3}\,{{\rm{Gpc}}}^{-3}\,{{\rm{yr}}}^{-1}\), consistent with ref. 40. The inferred \(\widetilde{m}\) is marked by the black symbol. Solid lines indicate the median merger rate and dashed lines the 10th–90th percentiles. A mass gap in the low-spin population, isotropic spins above \(\widetilde{m}\) and a sharp drop in the merger rate at the same mass value (the cliff) are all features consistent with a PISN gap populated by hierarchical mergers in dense star clusters.
The precise measurement of \(\widetilde{m}\) means that the mixture model is strongly favoured over models in which black holes of all masses share the same spin distribution. To quantify support for the transition-mass model, we compute the Bayes factor relative to the default model used by the LIGO–Virgo–KAGRA collaboration (see, for example, ref. 40), in which the full χeff distribution is described by a single truncated normal, \(P({\chi }_{\mathrm{eff}})={\mathcal{N}}({\chi }_{\mathrm{eff}};\,\mu ,\sigma )\), with no distinct high-mass component. We sample this hierarchical model using the same event set, likelihood and inference framework as in the main analysis. We obtain a Bayes factor B > 104 in favour of the model with a separate high-mass spin component.
The overall mass distribution shows several features. There are peaks at ~10 M⊙, 18 M⊙ and 38 M⊙, as also reported in ref. 40—these peaks were previously identified in ref. 41. There is a drop of nearly two orders of magnitude in the total merger rate at ~40 M⊙ (we name this feature ‘the cliff’ in Fig. 1). A rapid decline in the merger rate followed by a plateau (or a shallower decline) at \(\widetilde{m}\), as we found, is a generic feature of cluster formation models36,42. The plateau starts at the onset of PISNs and it is due to the emergence of binaries with components formed from previous mergers. According to these models, the edge of the transition between the two populations can be estimated from the primary-mass distribution alone as the value of maximum curvature above 30 M⊙ in the differential merger rate, Γ(m1) (‘The pair-instability mass gap in star-cluster models’). We measure this transition from the data at ~42 M⊙, consistent with the independent estimate based on the spin transition. A comprehensive interpretation of the cliff may require accounting for multiple concurring physical processes, as discussed in Supplementary Information.
The non-parametric analysis over a large dataset gives us confidence that a transition to a broader and more uniform distribution exists in the population. Motivated by this, we introduce a more informed parametric model in which the high-mass population is described by a uniform distribution with independent bounds, and use this model in what follows (unless otherwise specified). As before, the distribution below \(\widetilde{m}\) is well described by a normal distribution with positive mean, \(\mu =0.0{5}_{-0.02}^{+0.02}\), and small dispersion, \({\log }_{10}\sigma =-1.2{9}_{-0.17}^{+0.18}\). The distribution above \(\widetilde{m}\) is a broad distribution for which we infer upper and lower bounds of \({\chi }_{{\rm{eff,max}}}=0.4{9}_{-0.12}^{+0.14}\) and \({\chi }_{{\rm{eff,min}}}=-0.4{1}_{-0.35}^{+0.40}\), respectively. The lower bound is less well constrained than the upper bound, primarily due to selection effects and parameter-estimation uncertainties. Nevertheless, we find that χeff,min < 0 with 98.4% credibility, providing strong evidence for misaligned spins in the population. The median value \(\langle {\chi }_{{\rm{eff}}}\rangle =0.0{2}_{-0.17}^{+0.18}\) indicates that the distribution is statistically consistent with being symmetric around zero. Figure 2 shows the corresponding recovered χeff distributions.
Fig. 2: The χeff distributions below and within the PISN mass gap.
The alternative text for this image may have been generated using AI.
Thick black lines show the χeff distribution for the high-mass population, \(m > \widetilde{m}\), modelled as a uniform distribution with independent bounds, while thick red lines correspond to the χeff distribution of the low-mass population, \(m < \widetilde{m}\). Solid lines are the medians, while dashed lines show 10% and 90% of the distributions. The distribution for χeff in the PISN mass gap predicted under a hierarchical formation scenario is shown in green42. Light colour traces correspond to a single draw from the posterior distribution, providing a visual representation of the sample support from which the confidence intervals are constructed.
From the same parametric model we infer a characteristic mass scale of \(\widetilde{m}=44.{3}_{-3.5}^{+5.9}\,{M}_{\odot }\), a value close to the lower edge of the pair-instability mass gap reported in theoretical and numerical studies13,16,17,18. The posterior distribution of \(\widetilde{m}\) is shown in Fig. 3a.
Fig. 3: Constraints on the PISN transition mass and the 12C(α, γ)16O reaction rate.
The alternative text for this image may have been generated using AI.
a, Posterior distribution of the primary-mass value separating the two black hole populations, \(\widetilde{m}\). b, The corresponding posterior of S300. The latter is derived from \(\widetilde{m}\), using this as the value of the lower edge of the PISN mass gap.
Our results are consistent with a depletion of first-generation, low-spin black holes above \(\widetilde{m}\simeq 45\,{M}_{\odot }\). We reported a similar transition at \(4{6}_{-6}^{+7}\,{M}_{\odot }\) in GWTC-3 using 69 sources, with 11 having 90% of the m1 posterior distribution above 45 M⊙ after population reweighting36, and a transition to a higher-spin population was identified by others43,44. The consistent recovery in the new, larger catalogue containing 153 sources, with 34 above 45 M⊙, demonstrates that the feature is strongly driven by the data and it is not a statistical fluctuation. Constraints on the χeff distribution of the high-mass population are now tighter and more consistent with our theoretical interpretation. Using GWTC-3 we inferred \({\chi }_{{\rm{eff,max}}}=0.5{7}_{-0.19}^{+0.21}\) (uniform model), and could not decisively exclude χeff,max = 1 or χeff,max < 0 (ref. 45). With the expanded dataset we can instead rule out both χeff,max = 1 and χeff,max < 0 and infer χeff,max ≈ 0.5. This is notable because hierarchical mergers are the only astrophysical pathway that robustly predicts this upper bound36. The low-χeff tail is also better constrained. In GWTC-3, χeff,min and the left-hand tail were sensitive to prior choices (see, for example, Fig. 5 in ref. 45); with the larger sample we find −1 < χeff,min < 0 with substantially higher confidence, and a peak near the expected value χeff,min ≈ −0.5.
Together, these robust features indicate that the data are consistent with a hierarchical origin of the high-mass population: a mass gap in the low-mass/low-spin population, the onset of an isotropic and highly spinning population above \(\widetilde{m}\), the sharply defined upper bound of the χeff distribution at ~0.5 and the steep decline in the total merger rate (the cliff) near the transition. The transition mass itself, \(\widetilde{m}\simeq 45\,{M}_{\odot }\), matches stellar-evolution predictions for the onset of the pair-instability mass gap. Recent studies46,47 have identified a sharp decline in the merger rate—or possibly even a ‘gap’, as suggested by concurrent work48—for systems with secondary masses above ~45 M⊙, which is expected due to the rarity of binaries in which both components are second-generation black holes. Although PISN-gap mergers are expected mainly to involve a first-generation black hole and a merger remnant, we note that binaries where both black holes are merger remnants should also occur, as has been suggested for GW190521 (see, for example, ref. 49). Together, these independent lines of evidence favour hierarchical mergers as the most likely origin of the high-mass population. Such signatures are difficult to explain through isolated binary evolution, but arise naturally if the high-mass population is built from hierarchical mergers in dense stellar environments. With this level of confidence, we can now use this result to place direct constraints on massive-star evolution and the physics of the pair-instability process.
The location of the PISN boundary is ultimately set by stellar evolution physics, and in particular by the relative abundances of carbon and oxygen in the cores of very massive stars before collapse. These abundances depend on the 12C(α, γ)16O reaction rate, which governs the conversion of carbon into oxygen during helium burning13,16,17,18. A higher rate enhances oxygen production, leading to larger oxygen-rich cores and, consequently, to PISNs occurring at lower stellar masses. Conversely, a lower rate leaves behind more carbon, shifting the onset of pair instability to higher progenitor masses. Thus, measurements of the PISN mass gap from gravitational-wave observations of black hole mergers can provide an astrophysical constraint on the 12C(α, γ)16O cross-section, a quantity that remains one of the most important nuclear-physics uncertainties in massive stellar modelling15,50.
We assume that \(\widetilde{m}\) is the lower edge of the PISN mass gap, and follow refs. 15,51 to translate our inferred \(\widetilde{m}\) posterior into an estimate of the corresponding astrophysical S-factor at 300 keV, S300. The astrophysical S-factor rewrites a nuclear reaction cross-section by factoring out the strong Coulomb-barrier dependence, \(\sigma (E)=\frac{S(E)}{E}\,{{\rm{e}}}^{-2{\rm{\pi }}\eta }\), where E is the centre-of-mass energy and η the Sommerfeld parameter. We obtain \({S}_{300}=26{8}_{-116}^{+195}\,{\rm{keV\; b}}\) (90% credibility); we plot this probability distribution in Fig. 3. This estimate is consistent, within uncertainties, with recent nuclear-physics determinations52,53,54.
Our inference of the 12C(α, γ)16O S-factor from gravitational-wave data provides a novel, astrophysical constraint on a parameter that has long been central to stellar evolution theory. It relies solely on the assumption that the population with mass of ≳45 M⊙ consists entirely of second- (or higher-) generation black holes. Although direct nuclear-physics experiments have yielded estimates with large uncertainties15, our measurement achieves substantially tighter bounds, enabled by the sensitivity of the black hole mass spectrum to the details of helium burning. This improvement has wide-ranging implications: the carbon-to-oxygen ratio set by this reaction influences the core structure of massive stars, and thus affects the predicted rate of core-collapse supernovae, the maximum masses of neutron stars and the fate of red supergiants. It also governs the composition of white dwarfs (see, for example, ref. 55), with consequences for type Ia supernova explosions (see, for example, ref. 56), and shapes the nucleosynthetic yields that feed into galactic chemical evolution (see, for example, ref. 57). More broadly, the balance between carbon- and oxygen-rich material determines the conditions for planet formation and the likelihood of forming C-rich versus O-rich planetary systems (see, for example, ref. 58). Gravitational-wave astronomy therefore not only constrains the physics of compact objects, but also offers a new window into the nuclear processes that regulate stellar evolution and the chemical enrichment of the Universe.
We now use a non-parametric approach to search for additional isotropically spinning components in the data and to further test our result of a transition in spin properties at m1 ≃ 45 M⊙. We model the mixture fraction between the low-spin and the high-spin and isotropic populations as a non-parametric function of the primary mass. The posterior distribution of the mixture fraction is shown in Fig. 4, indicating that the fraction of isotropically spinning binaries is consistent with a sharp increase above ≳45 M⊙.
Fig. 4: Mass-dependent mixture fraction between the high- and low-spin populations.
The alternative text for this image may have been generated using AI.
a, The mixture fraction between the two black hole populations, modelled non-parametrically as a function of mass. The median and the 10% and 90% percentiles are shown. Light black lines are individual traces. b, The total differential merger rate as a function of primary black hole mass for reference. A transition above ~50 M⊙ is clearly recovered by this analysis. We also find a possible indication of an isotropic population at a primary mass of ~14 M⊙, although it is not statistically required by the data. This aligns with the mass dip in ref. 59 (green) and with GW241011 and GW241110, two O4b events interpreted as hierarchical mergers82 (pink). Together, these support our interpretation of a low-mass ‘valley’, which is populated by hierarchical mergers. Below ~10 M⊙, the posterior broadens, reverting to the prior due to the absence of sources there.
A new feature appears at m1 ≃ 14 M⊙, where the 90% bound of the mixing fraction rises to ~0.6. This coincides with a possible dip in the merger rate—previously identified in refs. 41,59. We interpret this as marginal evidence for an additional lower-mass gap in first-generation black holes that is populated by black holes formed from a previous merger60. While consistent with current data, this feature is not statistically required. Applying the same model to GWTC-34 yields an upper bound of ~0.16, showing that this feature only becomes discernible with the larger GWTC-4 catalogue.
The data indicate that nearly all primary black holes above 45 M⊙ involved in binary mergers possess high, isotropic spins. Explaining this within stellar evolution would require a mechanism that produces mass-dependent black hole spins at the end of massive-star lifetimes while also overcoming the pair-instability gap. The latter might be achieved through reduced stellar winds at low metallicity combined with the collapse of the residual hydrogen-rich envelope during a failed supernova26,61, but no explanation currently exists for the former. Another possibility is that stellar evolution could generate rapidly rotating black holes above ~45 M⊙ through fallback of angular-momentum-rich envelopes.
Our preferred explanation is that primary black holes above 45 M⊙ are the products of repeated mergers in globular clusters42. The inferred merger rate above this mass therefore provides a strong constraint on the initial cluster density: if all mergers above 45 M⊙ have this origin, the models of ref. 42 imply formation densities of ≳104 M⊙ pc−3. Repeated mergers may also occur in active galactic nucleus disks62 or nuclear star clusters32. Because these environments have higher escape velocities than do globular clusters, the detailed shape of the m1 distribution could help determine their relative contributions. After our work appeared on arXiv, ref. 63 likewise found a zero-symmetric χeff distribution above 45 M⊙. They showed that the two populations can also be distinguished by their mass-ratio distributions, with the higher-mass population favouring more asymmetric binaries, consistent with hierarchical formation in clusters. They noted, however, that the mass-ratio distribution of the high-mass population was not reproduced by their particular globular-cluster models. This could reflect either uncertainties in the cluster models and in the limited parameter space explored, or additional channels contributing within the PISN gap.
As the catalogue of detected binary black holes continues to expand with future observing runs, constraints on the pair-instability mass gap will sharpen, enabling increasingly stringent bounds on the 12C(α, γ)16O cross section. In the coming years, gravitational-wave population inference will thus not only elucidate the astrophysical environments where black holes form and merge, but also offer a new avenue to constrain fundamental nuclear reaction rates that underpin the evolution and fate of massive stars. At the same time, the identification of a population formed in dense star clusters offers a powerful opportunity to probe their initial conditions and evolutionary pathways across cosmic time.