Suspicious Siblings: The Distribution of Mass and Spin across Component Black Holes in Isolated Binary Evolution

Prior to the direct observation of compact binary coelescences via gravitational waves (GWs), the expected birth properties of black holes (BHs) relied primarily on highly uncertain predictions from stellar population modeling and a limited number of indirect observations. The past half-decade has brought an observational sample of binary black hole (BBH) mergers that has not only provided invaluable insights into the birth properties of BHs but also unveiled a number of unexpected systems that are in tension with the theoretical expectations of BH formation channels. Trends and features of the population as a whole are also becoming apparent, such as structure beyond a sharp cutoff in the BH mass distribution (e.g., Abbott et al. 2021a) and potential correlations between intrinsic parameters of BBH systems (e.g., Callister et al. 2021).

Certain BBH systems observed by the LIGO–Virgo–KAGRA (LVK) collaboration show unexpected features not just in their masses or spins individually, but their distribution of mass and spin across the two binary components. For example, GW190412 was the first BBH to have been measured with definitively asymmetric component masses (q ≃ 0.28, where q=m₂/m₁ with m₂ ≤ m₁) and nonzero spin (χ_eff > 0.14 at 90% credibility, where χ_eff=(m₁ a_1z + m₂ a_2z )/(m₁ + m₂) and a_z is the projection of the BH dimensionless spin aligned with the angular momentum (AM); Abbott et al. 2020). Because of the mass asymmetry in the signal, the component spin of the more massive primary BH was able to be disentangled from the leading-order spin term χ_eff, with the more massive BH having a dimensionless spin of a₁ > 0.22 at 95% credibility. More BBH systems with spinning component BHs have been identified in recent GW catalogs (e.g., Abbott et al. 2021b, 2021d, 2021e; Nitz et al. 2021; Olsen et al. 2022).

The combination of BH masses and spins across both binary components holds important clues about both the formation environment of BBH systems as well as the intricate physical processes occurring within and between their progenitor stars. In massive stars, AM transport between the stellar core and envelope significantly impacts the expected spin of the resultant BH. Models of AM transport in massive stars, such as the Taylor–Spruit magnetic dynamo (Spruit 2002), indicate AM transport to be highly efficient (Heger et al. 2005; Fuller et al. 2019). Under this paradigm, AM would thus be transported to the outer layers of the star during its giant phase and lost due to wind mass loss or Roche-lobe overflow (RLO) mass transfer (MT) onto a binary companion, resulting in a slowly spinning core and upon collapse a BH with a dimensionless spin extremely close to zero (Qin et al. 2018; Fuller & Ma 2019; though see Belczynski et al. 2020 for variations in this mechanism that can lead to slightly larger spins of a ≃ 0.1). The low effective spins in most BBHs also observationally hint at efficient AM transport in their progenitors (Abbott et al. 2021a, 2021c; Miller et al. 2020; Zevin et al. 2021), though certain systems are beginning to challenge the universality of quasi-isolated BHs having low spins at birth (e.g., Zevin et al. 2020a; Qin et al. 2022).

Even in the efficient AM transport paradigm, binary processes such as tides can induce high spins on stellar cores that can be preserved in their remnants. This can be accomplished in three general ways: (a) chemically homogeneous evolution, where two massive stars with near-equal masses in a tight binary system at zero-age main sequence (ZAMS) tidally interact on the main sequence, which induces strong rotational mixing that prevents expansion and AM loss, resulting in two massive, spinning BHs (De Mink & Mandel 2016; Mandel & De Mink 2016; Marchant et al. 2016); (b) Case A MT (i.e., while the donor is on the main sequence) for binaries on tight orbits of about a day or less, where the donor and accretor are tidally locked during MT and the envelope of the donor is stripped, leading the donor to never expand during its evolution into a Wolf–Rayet (WR) star (Qin et al. 2019); (c) tidal spin-up of a WR star (i.e., a naked helium star), either by an already-formed compact object following a stable or unstable MT event that hardens the binary to subday orbital periods (Detmers et al. 2008; Kushnir et al. 2017; Hotokezaka & Piran 2017; Zaldarriaga et al. 2018; Qin et al. 2018; Bavera et al. 2020; Olejak & Belczynski 2021; Steinle & Kesden 2021; Fuller & Lu 2022) or following a double-core common envelope (CE) scenario where two helium cores of supergiant stars are enveloped in the envelope of one or both of the supergiants and hardened through a CE phase to a point where the two cores can tidally interact (Brown 1995; Dewi et al. 2006; Hotokezaka & Piran 2017; Neijssel et al. 2019; Olejak & Belczynski 2021). Extremely metal-poor stars born in the early universe may also be able to retain a substantial hydrogen envelope and collapse into high-spinning BHs, though their contribution to the local merger population is likely small (Cruz-Osorio et al. 2021; Tanikawa et al. 2022). Following formation, a BH can also gain AM through accretion, though any appreciable spin-up would require the accretion to either be highly super-Eddington or transpire over timescales far longer than the evolutionary timescales of massive stars (van Son et al. 2020; Bavera et al. 2021b; Qin et al. 2022).

Systems with unequal masses and spinning primaries provide a challenge to the isolated binary evolution scenario. BBHs that result from chemically homogeneous evolution strongly favor near-equal-mass systems (Marchant et al. 2016; Mandel & De Mink 2016; du Buisson et al. 2020). The Case A MT scenario has been invoked for explaining the large inferred spins of BHs in high-mass X-ray binaries (Qin et al. 2019), though binary modeling finds that the high-mass X-ray binaries in the Milky Way are unlikely to form merging compact binary systems (Belczynski et al. 2011, 2012; Neijssel et al. 2021). The BH–WR tidal spin-up scenario is predicted to be common for post-CE binaries (Bavera et al. 2021b), though the efficiency of this pathway is highly dependent on CE ejection efficiency and can only impart spin on the second-born BH progenitor. Though the tidal spin-up of two WR stars following a double-core CE may lead to significant spins in both BHs (Olejak & Belczynski 2021), it typically leads to near-equal-mass mergers and operates at a much lower rate than the BH–WR tidal spin-up scenario (Neijssel et al. 2019). Lastly, although spinning up the first-born BH through accretion can be accomplished via highly super-Eddington MT, if accretion is pushed too high, the rate of mergers from this channel can drop by orders of magnitude (Bavera et al. 2021b).

The BH–WR tidal spin-up scenario described above may be a valid explanation for the observed spins of primary BHs if the second-born BH can oftentimes be more massive than the first-born BH. In isolation, a more massive star at ZAMS will typically lead to a more massive remnant (though see, e.g., Patton et al. 2022). However, binary interactions throughout the coevolution of the two stars can alter this picture. In particular, since the more massive star at ZAMS will almost always evolve off the main sequence first, it will be the first to overflow its Roche lobe and transfer material onto its companion. The fraction of transferred mass deposited onto the companion depends on an uncertain MT accretion efficiency (Bouffanais et al. 2021b), and can potentially lead to a mass ratio reversal (MRR), where the star that was originally less massive at ZAMS has its mass inflated from the accreted material and leads to a more massive remnant (e.g., Olejak & Belczynski 2021; Hu et al. 2022). The originally more massive star will still proceed through the remainder of its evolution quicker, and create the first compact object in the binary. Then, the originally less massive star can be stripped of its envelope and harden its orbit with the compact object companion during the second MT episode, and be tidally spun up as a WR star.

In this paper, we explore the viability of the MRR/tidal interaction and BH accretion spin-up scenarios for forming the asymmetric-mass, spinning-primary BBH systems observed by the LVK. In Section 2, we overview our population models and the determination of BH spins through tidal spin-up and accretion. Analysis of these models, with a particular focus on mass ratios, MRR, and BH spins, is in Section 3. We contextualize our results with GW events and population properties in Section 4, and discuss the broader implications and caveats in Section 5. Throughout this work, we assume solar metallicity of Z_⊙ = 0.017 (Grevesse & Sauval 1998) and Planck 2018 cosmological parameters (Alves et al. 2020).

We use the open-source binary population synthesis code COSMIC ⁵ (Breivik et al. 2020) for modeling the populations of BBH mergers. COSMIC is based on single-star evolutionary tracks from the SSE code (Hurley et al. 2000) and the binary star implementations of BSE (Hurley et al. 2002). A large number of updates have been made to the physical prescriptions used for initial conditions, winds, the onset of MT and system evolution during RLO, remnant formation, and natal kicks. Rather than simulating a fixed number of systems with specific physical assumptions at a given metallicity, COSMIC samples the systems until user-specified convergence criteria have been reached in the population. ⁶ Although we highlight physical assumptions pertinent to this work in the following section, we refer the reader to Breivik et al. (2020) for further details.

2.1. Physical Assumptions

A large number of physical uncertainties embed models of massive-star binary evolution, which have significant impacts on the population properties of their BBH remnants (e.g., Giacobbo et al. 2018; Giacobbo & Mapelli 2018; Kruckow et al. 2018; Bouffanais et al. 2019; Stevenson et al. 2019; van Son et al. 2020; Bavera et al. 2021b; Belczynski et al. 2022; Santoliquido et al. 2021; Zevin et al. 2021). Although recent studies have become increasingly thorough in their coverage of this highly uncertain parameter space (see Broekgaarden et al. 2021 for a recent overview), it is computationally expensive and can make interpretation difficult. We instead narrow our parameter space coverage to uncertainties that have the strongest effect on MRR and BH spin-up.

The main parameter that is varied between models is the accretion efficiency:

where is the mass-loss rate of the donor during RLO and is the mass accepted by the accretor. The additional mass is lost from the system with an AM as if it were a wind from the accretor. We simulate five variations for f_acc: [0.0, 0.25, 0.5, 0.75, 1.0] where f_acc = 0.0 means that the accretor accepts no mass from the donor, and f_acc = 1.0 means that the accretor accepts all of the mass transferred by the donor. Although the timescales pertinent to the response of the envelope to MT may provide a more realistic description of the amount of material a star can accrete in a given amount of time, this simple parameterization can capture the extreme limits of accretion efficiency, which has an important impact on possible MRR during stable MT and subsequent tidal spin-up of the second-born BH. This parameter also impacts the compact binary merger rates and mass spectra (see, e.g., Kruckow et al. 2018), and has the promise of being constrained with GW data (Bouffanais et al. 2021b). The mass-loss rate from the donor during RLO, , follows the prescription of Hurley et al. (2002).

In addition, we simulate two separate assumptions for the amount at which BHs can accrete material above the Eddington rate, which impacts the accretion-induced spin-up of the first-born BH:

where R_s is the Schwarzschild radius of the BH in units of solar radii and γ_Edd is a multiplicative factor that allows for super-Eddington accretion (i.e., γ_Edd = 1 would limit accretion onto a BH to the Eddington rate). The maximum rate at which a BH can accrete is thus given by

Note that we allow f_acc to alter the accretion efficiency of compact objects as well as (nondegenerate) stars in our models, and thus at f_acc = 0, BHs will not gain mass through accretion. We choose two values for γ_Edd: [1, 10⁵]. The highly super-Eddington parameterization of γ_Edd = 10⁵ is extreme but chosen because it can lead to a significant number of BHs to be spun up through stable MT onto the BH, whereas lower values for γ_Edd are much less capable (see, e.g., Bavera et al. 2021b).

The final parameter variations we explore govern the evolution through a CE phase. These affect both the survival of systems that evolve through a CE phase and post-CE orbital separations, thereby impacting the ability of the progenitors of the second-born BH to tidally spin up. We model three variations for the efficiency at which orbital energy of the inspiraling binary is transmitted into the energy needed to eject the envelope, assuming the α − λ formalism for energetics during CE evolution (Webbink 1984): α_CE = [0.5, 1.0, 2.0]. We use the variable prescription in Claeys et al. (2014) for determining the value of λ. Stellar type-dependent mass ratios that determine whether MT proceeds stably or unstably follow Neijssel et al. (2019), though we also run a subset of models that follow Belczynski et al. (2008) to investigate the impact this parameterization has on MRR and tidal spin-up (see Appendix B).

All models sample the binary initial conditions independently following Kroupa (2001) and Sana et al. (2012), and assume a pessimistic CE scenario in which Hertzprung gap stars that experience unstable RLO merge in the CE and do not form compact binaries (see Belczynski et al. 2008). We assume a fixed binary fraction of f_bin = 0.7, which only affects the arbitrary normalization during the resampling described in Section 2.3. Wind mass loss follows Hurley et al. (2002) with updates for O and B stars (Vink et al. 2001) and metallicity-dependent WR winds (Vink & de Koter 2005). Remnant masses are determined assuming the delayed supernova engine prescription of Fryer et al. (2012) with neutrino mass loss as implemented in Zevin et al. (2020b). A maximum neutron star mass (and therefore minimum BH mass) of 3 M_⊙ is assumed. Supernova kicks, which can tilt the orbital plane and lead to BH spins that are misaligned from the orbital AM, are drawn from a Maxwellian with a dispersion of 265 km s⁻¹ (Hobbs et al. 2005) and are fallback modulated, leading to a suppression of kick strength as a function of BH mass (e.g., Rodriguez et al. 2016). BH masses as a result of pulsational pair instability and pair instability supernovae are treated using fits to the results of Marchant et al. (2019).

Population models only retain systems that result in a BBH that merges within a Hubble time. For each population model, we simulate 12 fixed metallicities spaced uniformly in log between Z = [10⁻⁴, 0.03]. The permutation of all these model variations leads to 5 [f_acc] × 2 [γ_Edd] × 3 [α_CE] × 12 [Z] = 360 individual population sequences. Throughout this work, we present the results that combine the 12 discrete metallicities for each model set according to the joint star formation and metallicity evolution of the universe (see Section 2.3).

2.2. BH Spins

Tracking the spin evolution of BH progenitors is difficult in rapid population synthesis, and typically relies on simplified prescriptions. This is because rapid population synthesis codes lack information about the internal structure of a star, and thus cannot accurately model the AM transport between binary components or back-reaction on the structure and evolution of each individual star. We assign spins to the first-born BH (a_1b) and second-born BH (a_2b) in our population during post-processing, assuming that AM transport is highly efficient (Spruit 2002; Fuller & Ma 2019) such that the first-born BH can only attain a nonzero spin through accretion after it becomes a compact object, and the second-born BH can only attain a nonzero spin through tidal spin-up of its WR progenitor. Thus, BHs born in quasi-isolation where their progenitors are not susceptible to tidal spin-up are assumed to be Schwarzschild BHs, though we examine the impact of relaxing this assumption on our results to account for other possible mechanisms for inducing spin in Section 5. We note that post-processing spins in this way leads to a slight inconsistency in our determination of inspiral times through GW emission, as spins aligned with the orbital AM will slightly expedite the eventual merger. However, this has a minor effect for the timescales considered, because spin effects enter at a higher post-Newtonian order and do not impact the inspiral rate significantly until the binary is close to merger (Damour 2001).

The spin of the first-born BH is calculated assuming stable MT from a disk of material being accreted from the innermost stable circular orbit of the BH as in Thorne (1974):

where M_i is the mass of the BH prior to MT and M_f is the mass of the BH after MT has ceased. For the determination of a_1b, we only consider the mass gained from stable RLO MT and neglect any potential mass gain from wind accretion. We also assume that the BHs gain no mass when inspiraling through a CE (though see, e.g., MacLeod & Ramirez-Ruiz 2015; Cruz-Osorio & Rezzolla 2020). The spin of the second-born BH is determined using the semianalytic fits from Bavera et al. (2021a), which are based on the detailed spin evolution of BH–WR systems during tidal spin-up using the MESA simulations (Paxton et al. 2011, 2013, 2015, 2018, 2019) under the POSYDON ⁷ framework (Fragos et al. 2022). Bavera et al. (2021a) found the spin of the second-born BH to be well approximated by a quadratic function dependent on BH–WR log-orbital period , which implicitly is dependent on the mass of the WR star M_WR at helium depletion:

where with coefficients , , and determined through least-square minimization; see Bavera et al. (2021a). Since the grid of systems used for the fit only extended down to orbital periods of 0.1 day (given by the physical limit of RLO at zero-age helium main sequence in the MESA simulations), for systems in our population with p < 0.1 day, we assume their orbital periods are 0.1 day when computing the fit. Across our population models, 0.2%–2.9% of BH–WR systems have orbital periods of p < 0.1 day and are extrapolated in this manner.

We show the second-born BH spins for a single population model in Figure 1. For M_WR < 14 M_⊙ at low periods (p < 0.3 day), the resultant BH spin decays as a function of WR star mass. Such a trend is dictated by the core-collapse mechanism. For stars with carbon-oxygen core masses , the Fryer et al. (2012) delayed prescription predicts the mass ejection in the supernova event. The ejected mass will carry away the AM stored in the outer layers, lowering the fraction of AM transferred to the BH from the WR star. In the simulations of Bavera et al. (2021a), we also see a similar BH spin decay for BH–WR systems with M_WR >40 M_⊙ and low orbital periods (p < 0.3 day). Such suppression is due to strong WR stellar winds and pulsations due to the pair instability process, which both cause the star to lose a nonnegligible fraction of material from its envelope, depleting the WR AM reservoir.

Zoom In Zoom Out Reset image size

One additional subchannel of BBH formation from isolated binary evolution that is present in our populations is the double-core CE channel, in which the helium cores of two supergiant stars proceed through a CE in one or both of the envelopes of the stars, leading to a tight WR–WR system where both binary components can be tidally spun up. Since this channel is subdominant and we currently lack the simulation grids and fits to monitor this subdominant tidal spin-up scenario, we exclude such systems from our population and focus solely on the standard CE and stable MT channels that lead to tidal spin-up. This channel also typically forms near-equal-mass BBH systems and therefore does not strongly affect our key results. We comment more on our exclusion of this subchannel and implications in Section 5.

2.3. Star Formation and Metallicity Evolution

As described in Section 2.1, for each model (parameterized by f_acc, γ_Edd, and α_CE), we simulate 12 discrete metallicities. Each binary is evolved for a full Hubble time t_H. To get a population of merging BBHs that is representative of the underlying population in the universe, we resample 2 × 10⁵ binaries from the metallicity runs and populate the binaries in redshift according to star formation history and metallicity evolution as predicted by the Illustris-TNG simulations (Nelson et al. 2015).

From Illustris-TNG, we have a grid of stellar mass formed over metallicity and lookback time in a 100 comoving Mpc³ cube. Marginalizing over the metallicity axis provides the star formation rate density as a function of redshift, ψ(z). However, although our population modeling tracks the total stellar mass sampled across all single stars and binaries, our population models only retain systems of interest (i.e., those that form BBHs). Therefore, drawing the birth redshifts for our target population according to this distribution does not account for the relative formation efficiency of our target population at each redshift, which is dependent on the metallicity distribution at each redshift. For each metallicity model Z_sim, we have a formation efficiency

where N_BBH is the number of BBHs formed and M_sample is the total mass sampled in the particular metallicity model, accounting for the entirety of the initial mass function and binary fraction. At a given redshift z, we can construct a distribution of metallicities from the cosmological model. The cumulative distribution function of this distribution can be used to determine the relative contribution of each discrete metallicity at a given redshift. We define p_Z (Z∣z) as the support of the metallicity cumulative distribution function at redshift z closest to a metallicity model Z_sim in log-space, such that . Combined, gives the number of systems per sampled stellar mass that should be drawn from each metallicity model at a given redshift.

The number of BBH systems formed per comoving volume at a particular discrete metallicity across all redshifts is thus

where we perform the change of variable with . This integral gives us the total number of BBH systems formed per source-frame year in a comoving box at a discrete metallicity over the entire history of the universe. At a given metallicity, we densely discretize redshifts and evaluate this integral between the bounds of our redshift bins, which gives a relative weight at a particular redshift and metallicity that we use when resampling our population models.

With the relative probability of drawing a system from our target population at a given redshift and metallicity in hand, we randomly draw the birth redshifts and metallicities of 2 × 10⁵ systems for each population model set. We remove the systems from our resampled population that merge after the present day; that is, systems that satisfy the condition τ(z) − t_delay < 0, where t_delay is the delay time (defined as the time between ZAMS and BBH merger) and τ(z) is the lookback time of a system born at redshift z. This eliminates 4%–13% of systems from our populations of merging BBHs. The codebase used for generating a metallicity and redshift resampled population based on the COSMIC models using various star formation history and metallicity evolution assumptions is available on Github. ⁸

We now investigate the viability of isolated evolution at forming systems with asymmetric masses and spinning primaries. We first look at the component spin distributions and mass ratio distributions individually for our array of population models to identify trends across our physical assumption variations in Sections 3.1 and 3.2. We then look at broader population properties in mass ratio–spin space in Section 3.3.

3.1. Spin Magnitudes of the First- and Second-born BH

Under the assumption of efficient AM transport in massive stars, spin in the first-born BH can be achieved through accretion via stable MT onto the already-formed BH by a stellar companion. Figure 2 shows the distribution of first-born BH spins a_1b across our various physical assumptions. The distributions transition from a flat region, which indicates the branching fraction of stable MT systems in each model, to a monotonically increasing behavior; the transition at a_1b ∼ 3 × 10⁻⁴ in all models corresponds to the minimum accretion timescale found in our populations of ∼0.01 Myr. Appreciable spin in the first-born BH can only be attained via accretion if the BH can accrete far above the Eddington limit. For Eddington-limited accretion onto the first-born BH, we find that no more than 0.9% of systems are spun up beyond a_1b > 0.01 across all population models. This is expected due to the short evolutionary timescales of the massive stars that are BH progenitors, which have post-main-sequence lifetimes of ≲1 Myr and thus can only accrete at most a fraction of a solar mass if accretion is limited to the Eddington rate. More systems proceed through this stable MT channel at lower accretion efficiencies (≃8% at f_acc = 0.25 compared to ≃2% at f_acc = 1 in the Eddington-limited models), though higher values for the accretion efficiency lead to larger typical first-born spins when only considering the stable MT channel. Variations in the CE efficiency, shown with different linestyles in Figure 2, have a minor and indirect impact on the first-born spin distributions; these slight variations are due to the changes in the relative fraction of merging BBHs that proceed through a CE phase rather than stable MT.

Zoom In Zoom Out Reset image size

Increasing the accretion limit onto BHs to 10⁵ times the Eddington rate drives the first-born BH spins to more extreme values. The differences in the transition between a flat behavior and monotonically increasing behavior for different accretion efficiency models in Figure 2 indicates that the maximum MT rate is now often imposed by the f_acc rather than the γ_Edd, with the transition happening at larger spins as f_acc increases. For the nonzero accretion efficiencies f_acc > 0, a sizeable fraction of systems in our populations have first-born spins with a_1b > 0.5: ≈ 6% in our f_acc = 0.5 models. Almost no systems have first-born spins of a_1b < 0.5 in any of our super-Eddington accretion models. Though we do not explicitly consider BBH merger rates here, we note that significantly increasing the possible accretion rate onto BHs drives down the expected merger rate of systems with highly spinning first-born BHs due to the conservative MT not shrinking the orbit as efficiently as nonconservative MT (Bavera et al. 2021b). This is seen in the bottom panel of Figure 2; though increasing the accretion efficiency does indeed lead to more highly spinning BHs from the stable MT channel, the relative contribution of this channel compared to the full BBH population decreases.

In Figure 3, we show the spin distributions for the second-born BHs a_2b in our models, which is driven by tidal spin-up of He cores by the first-born BH. For all models, we find a broad range of spin magnitudes ranging from nonspinning to maximally spinning, with little dependence on γ_Edd because the CE channel dominates over the stable MT channel in our models. The peak at a_2b ≃ 0 in all models is the result of systems that have orbital periods outside of the spin-up regime at WR–BH formation; see Figure 1. Across our models, we find that 22%–52% of secondary spin magnitudes are a_2b < 0.05.

Zoom In Zoom Out Reset image size

We find a general trend of second-born spin distributions pushing to larger values with increasing accretion efficiency and decreasing CE efficiency, in agreement with Bavera et al. (2021b). For f_acc = 0.5 and α_CE = 1, 31% of systems have second-born BH spins of a_2b > 0.5. This fraction drops to 18% for α_CE = 2 and increases to 52% for α_CE = 0.5. The increase in high-spinning BHs with decreasing CE efficiency is due to the less efficient CE phases leading to more hardened orbits post-CE, which makes the stripped He core more susceptible to tidal spin-up. At α_CE = 1, the fraction of systems with a_2b > 0.5 drops to 18% for f_acc = 0 whereas at f_acc = 1 it increases to 33%. This trend of increasingly high-spinning second-born BHs with increasing accretion efficiency is due to the progenitor of the second-born BH gaining more mass during the first stable MT phase; the post-CE separation scales inversely with the product of the total mass and envelope mass of the donor pre-CE.

3.2. Mass Ratios and Mass Ratio Reversal

We now turn to the mass ratio distributions in our populations. Since the order at which BHs are born in a binary determines the mechanism responsible for their potential spin, we choose to define the birth-mass-ratio parameter q_b = m_2b/m_1b, where m_1b and m_2b are the masses of the first-born BH and second-born BH, respectively. Thus, for values q_b > 1, the system proceeded through a MRR that led to the less massive star at ZAMS forming in the more massive BH in the binary. Note that in our models the only mechanisms for inducing the spin on BH remnants are accretion onto the BH and tidal spin-up, and thus the first-born BH is only susceptible to the former, and the second-born BH is only susceptible to the latter of these processes.

We show the population distributions of the mass ratio parameter q_b in Figure 4. The accretion efficiency f_acc has the largest impact in determining whether MRR ensues, as this parameters governs the amount of material the less massive star at ZAMS can accept from the initially more massive donor after it evolves off the main sequence and overflows its Roche lobe. For example, in our default model with α_CE = 1 and γ_Edd = 1, we find 0.3% of merging BBH systems to have q_b > 1 for f_acc = 0, which rises to 47% at f_acc = 0.5 and 72% at f_acc = 1. Despite this, the second-born BH rarely reaches masses of more than ≃3 times the mass of the first-born BH; in our most optimistic MRR models, we find the birth-mass-ratio parameter to be constrained to q_b < 2.2 at 99% confidence, and it never reaches a value larger than q_b = 3.7.

Zoom In Zoom Out Reset image size

Increasing the accretion limit onto BHs leads to more systems occupying the low-end tail of the q_b distribution. This is due to the first-born BH being able to gain significantly more mass from accretion relative to the Eddington-limited case, and can generate systems with mass ratios of q_b ≃ 0.09 in the most extreme cases. However, these extreme mass ratio systems are still rare in the context of the full population; at f_acc = 0.5 and α_CE = 1, we find that 99% of systems have q_b > 0.18 for γ_Edd = 10⁵, which raises to q_b > 0.33 for the Eddington-limited case with γ_Edd = 1. Thus, when considering all systems regardless if they proceeded through a MRR or not, we find that forming binaries with mass asymmetries of more than 3:1 is a rare occurrence unless super-Eddington accretion is invoked.

3.3. Population Properties

We now examine both component spins and mass ratios simultaneously to make broader statements about the viability of generating asymmetric-mass systems with spinning primaries. Figure 5 marks the fraction of systems in our full BBH population that satisfy particular criteria of interest. We specifically highlight the systems that proceed through MRR (q_b > 1), systems that have significant asymmetries in their BH masses (q_b > 2 or q_b < 0.5), and systems that achieve a significant primary BH spin (a > 0.2).

Zoom In Zoom Out Reset image size

As shown in Figure 4, the number of systems that proceed through a MRR is significantly enhanced with increasing accretion efficiency (blue and green lines). However, even for perfectly conservative MT where all of the material from the donor is accepted by the accretor, the number of systems that have significantly asymmetric masses and nonnegligible spins in the primary BH is small (yellow and orange lines). Across our Eddington-limited population models, 0.6%–1.5% satisfy this criteria at f_acc = 1; this number drops by an order of magnitude or more for lower values of f_acc. Therefore, in the efficient AM transport paradigm, we find the formation of systems with asymmetric masses and spinning primary BHs to be inefficient even for the most optimistic assumptions about MT physics.

For our models that allow for highly super-Eddington accretion (γ_Edd = 10⁵), the percentage of systems with nonnegligible primary BH spins and mass ratios greater than 2:1 can exceed 8% and drop to less than 0.5% depending on the choice of accretion efficiency, so long as f_acc > 0. When the accretion onto BHs is Eddington limited, no systems in our population meet this criteria since the first-born BH spin-up is highly suppressed (see Figure 2). Thus, in either spin-up scenario, we find the fraction of systems that have significant mass asymmetries (>2:1) and significant primary spins (>0.2) to be percent-level at best, and requires highly conservative MT or highly super-Eddington accretion to achieve these numbers.

In the previous section, we examined the efficiency of forming systems with asymmetric masses and spinning primary components, and found that such systems are relatively rare occurrences in the underlying population even when taking optimistic assumptions about MT efficiency. We now turn to the implications of these results, in the context of explaining both the individual events and general trends that are observed in the BBH population.

A number of systems in the current population of BBHs exhibit interesting combinations of mass ratio and spin. For example, GW190412 and GW190517_055101 are both high-confidence events with positive effective spins of and , respectively, at 90% credibility (Abbott et al. 2021b). The primary BH in GW190412 is 3–4 time more massive than the secondary, with a mass ratio of (Abbott et al. 2020), whereas the mass ratio of GW190517_055101 is more consistent with unity (, Abbott et al. 2021b). Because of the large measured effective spins and/or mass asymmetry, the spin magnitudes of the primary BH can be constrained to and for GW190412 and GW190517_055101, respectively, whereas the spin of the secondary BH is unconstrained. ⁹ The extended catalog for the first half of the third observing run uncovered the system GW190403_051519, and although it has a relatively low astrophysical probability compared to GW190412 and GW190517_055101 and a primary mass that is in tension with the maximum BH mass expected from the pair instability process, its significantly asymmetric masses () and high primary spin () make it a system of interest. In the most recent GW catalog published by the LVK, the high-confidence event GW191109_010717 also has intriguing spin signatures with the primary BH spinning at (Abbott et al. 2021b, 2021e), although its high mass and in-plane spin may be indicative of a dynamical formation scenario. Other GW catalogs outside the LVK have also uncovered additional systems with interesting configurations of spin and mass ratio that may help to probe the viable spin-up mechanisms (e.g., Nitz et al. 2021; Olsen et al. 2022).

These BBHs, particularly GW190412, are potential exemplars of systems whose progenitors underwent a MRR where the more massive BH was born second and its progenitor was tidally spun up (Abbott et al. 2020; Zevin et al. 2020a; Olejak et al. 2020; Olejak & Belczynski 2021), or alternatively having the first-born primary BH spun up by super-Eddington accretion during the second stable MT episode (Bavera et al. 2021b). As we showed in Figure 3, population modeling predicts a significant number of systems that can have appreciable spins across the entire range of physical spin magnitudes in the underlying population. However, the formation of systems that can also achieve spinning primary components drops precipitously as the mass of the primary relative to the mass of the secondary increases, as shown in Figure 5. Although the mass ratio distribution of GW190517_055101 is too broad to definitively rule out this formation mechanism, we find the formation of GW190412 through MRR and tidal spin-up to be improbable. Taking the upper and lower values of GW190412's 90% symmetric credible interval for mass ratio and primary spin, respectively (i.e., taking values closest to nonspinning and equal mass), we find that <0.1% of systems satisfy this criteria in our Eddington-limited accretion model even when we take the optimistic assumption of f_acc = 1. For our highly super-Eddington models, we find the formation of such systems to be more likely, with up to 8% of systems satisfying this criteria.

We find the detection of asymmetric-mass, spinning primary BH systems through the tidal spin-up scenario to be less likely when selection effects are accounted for, consistent with the findings of Bavera et al. (2021b). Since the tidal spin-up of a WR star by a BH requires subday orbital periods (see Figure 1), these systems have short inspiral times when they form BBHs. Systems will more readily proceed through tidal spin-up at lower metallicities (and therefore higher redshifts) because the WR stellar winds are weaker, which lessens orbital expansion and keeps the BH–WR system in the regime where tidal spin-up can still be efficient (Bavera et al. 2020). Therefore, many tidal spin-up systems will merge outside of the horizon of current ground-based GW detectors, and a lower fraction of systems in the detectable population is predicted to have high spins from tidal spin-up relative to the underlying population (e.g., Safarzadeh & Hotokezaka 2020; Bavera et al. 2021b). Applying selection effects to our population causes the probability of detecting systems similar to GW190412 through the MRR and tidal spin-up channel to further decrease by more than an order of magnitude. However, we find the super-Eddington accretion spin-up scenario to be less impacted by selection effects, still accounting for a percent-level number of systems in the detectable population. We describe our implementation of selection effects in Appendix A.

In addition to the individual BBH systems in the GW-detected population having interesting parameter estimates, populations analyses are useful for uncovering trends and features in the data as a whole (see, e.g., Mandel et al. 2019; Vitale et al. 2020). One example that is particularly pertinent to mass ratios and spin configurations is that the population of BBHs shows tantalizing signs of (anti)correlations between mass ratio and effective spin (Callister et al. 2021; Abbott et al. 2021c). In Figures 6 and 7, we show joint distributions of the birth-mass-ratio parameter, q_b, and effective spin, χ_eff, in our BBH populations for γ_Edd = 1 and γ_Edd = 10⁵, respectively. A few GW events with interesting mass ratio and spin measurements are overplotted for comparison.

Zoom In Zoom Out Reset image size

We find no strong indication of trends in effective spins as a function of mass ratio (for either q_b > 1 or q_b < 1) in our Eddington-limited population models shown in Figure 6. The effective spin distributions tend to narrow closer to χ_eff ≃ 0 as mass ratios decrease from q_b = 1 to q_b = 0. For values of q_b > 1, the effective spin distributions remain broad, but retain a peak in the distribution at χ_eff ≃ 0.

This situation differs in the super-Eddington models, as shown in Figure 7. Significantly increasing the accretion limit onto BHs leads to the stable MT channel populating the region of parameter space with large effective spins and small mass ratios (q_b ≪ 1). The combination of the bulk of systems dominating the low-spin and equal-mass regime, MRR and tidal spin-up accounting for mildly asymmetric masses and moderate spins, and stable MT with super-Eddington accretion producing low mass ratios and high spins can lead to an apparent trend between mass ratio and effective spin, and a potential means of generating the anticorrelation observed in GW data. Visually, the super-Eddington stable MT channel can populate the regions of parameter space occupied by unique systems such as GW190412. We also see an increased number of systems with q ≪ 1 and antialigned spins, which result from the second-born BH being lower mass and more susceptible to strong natal kicks that can significantly tilt the orbital plane. Although the relative fraction of BBHs formed via CE and stable MT channels may not be accurately captured by rapid population synthesis approaches (see Section 5 for further discussion), the combination of these isolated evolution subchannels may provide one explanation for anticorrelations between mass ratio and effective spin. However, our models find that extreme assumptions regarding the accretion rate onto BHs are required for such an interpretation.

Zoom In Zoom Out Reset image size

It is interesting to note that certain analysis techniques can have a significant impact on measurements of spins in individual BH systems and the population of component spins. For example, rather than inferring the properties of the most massive and least massive BHs in a particular system, Biscoveanu et al. (2021) show that, by instead inferring the parameters of the highest-spinning and lowest-spinning components of the system, spin magnitude measurements can be more informative especially for near-equal-mass systems. Such variations to the standard analyses may be useful in better constraining the component spin magnitudes, though it yields less differences in parameter estimates for unequal-mass systems such as those considered in this work, due to the spin of the heavier BH playing a larger role in measured spin quantities such as the effective spin. Nonetheless, such techniques may uncover more component BHs with significant spin magnitudes in the broader population (Abbott et al. 2021c).

In this work, we simulated a suite of population models to investigate the joint distribution of masses and spins across an astrophysical population of BBHs. In particular, we focus on the ability of stellar progenitors to proceed through a MRR where the lighter star at ZAMS leads to the more massive compact object, and the potential means of generating spin in BHs even with efficient AM transport in their stellar progenitors. Our key results are as follows.

• 1.  
  Although the distribution is peaked near zero, the second-born BH can attain a range of spin magnitudes from near zero to maximally spinning through tidal spin-up, with less efficient CE phases leading to larger typical spins. Across our model variations, 4%–57% of systems have second-born BHs with spin magnitudes >0.5 (Figure 3).
• 2.  
  Assuming highly efficient AM transport and Eddington-limited accretion, the spin of the first-born BH will naturally be small, with typical values ≪10⁻². However, increasing the accretion rate to many orders of magnitude above the Eddington limit will lead to a sizeable fraction of first-born BHs with significant spin. Though higher accretion efficiencies lead to larger typical spins from this subchannel, it causes the contribution of this subchannel to drop relative to the full population (Figure 2).
• 3.  
  Naturally, the number of systems that proceed through MRR is highly sensitive to the assumed accretion efficiency at which the accreting star accepts material from the donor star. Our parameterization choices for the accretion efficiency bracket the range of systems that have a more massive second-born BH to 0.3%–72%, though even at a low accretion efficiency of f_acc = 0.25, more than 18% of systems have a more massive second-born BH.
• 4.  
  It is rare for systems to form with asymmetric masses and spinning primary BHs through MRR and tidal spin-up; across our various model assumptions with f_acc > 0, 10⁻⁵ − 10⁻² of systems in the underlying population have second-born BHs that are twice as massive as their first-born counterparts and have spin magnitudes >0.2. These numbers decrease by an order of magnitude or more when GW-detector selection effects are considered due to the short delay times of this channel and lower efficiency of tidal spin-up at lower redshifts and higher metallicities. The number of systems where the first-born BH is twice as massive as the second-born BH and spinning at >0.2 is negligible for our Eddington-limited models, though can account for 5 × 10⁻³ − 8 × 10⁻² of systems in the underlying population of our models that assume highly super-Eddington accretion. We therefore find that the formation of BBH systems with asymmetric masses and spinning primary BHs, such as GW190412, is unlikely through the channels considered in this work unless highly super-Eddington accretion is invoked.
• 5.  
  We find no indication of a correlation between mass ratio and effective spin in our model variations other than in the models that assume highly super-Eddington accretion onto BHs.

Though we aim to survey the parameterizations of binary stellar evolution that most impact MRR and tidal spin-up, the vast array of physical uncertainties that embed binary population modeling prevent us from being completely exhaustive in our coverage of parameter space. We anticipate that the criteria determining the onset of unstable MT could lead to potential changes in our results. This is typically encoded in rapid population modeling via an array of stellar-type-specific critical mass ratio values, , where MT will be unstable if the donor mass is sufficiently more massive than the accretor mass. One variation in this parameter for a subset of our other model variations is described in Appendix B. However, it is possible that the values typically used in population synthesis are not representative of the true physical picture. For example, Gallegos-Garcia et al. (2021) found that rapid population synthesis may severely overpredict the number of systems that proceed through a successful CE phase and form BBH systems. In Figure 7, the high-density region at low mass ratios (q_b ≪ 1) and large, positive effective spins is populated almost exclusively through systems that do not proceed through CE evolution, and could potentially lead to a mass ratio–effective spin correlation in the full BBH population if the CE channel is suppressed, though this feature is only apparent in the models that assume super-Eddington accretion onto BHs (see also Bavera et al. 2021b). In our Eddington-limited models, we find that the stable MT channel typically does not harden BH–WR binaries to orbital periods of less than 1 day, and thus is inefficient at spinning up the second-born BH progenitor through tides. However, we note that the work investigating tidal spin-up in isolated evolution with other population synthesis codes, such as Olejak & Belczynski (2021), finds that evolutionary sequences that do not involve a CE phase can still lead to an appreciable fraction of systems that are spun up via tides.

One subchannel present in our models, though excluded in our analysis, is the double-core CE channel, in which both helium cores can tidally interact and spin up before collapsing into BHs. This channel is excluded (a) because we lack detailed simulations grids following tidal spin-up in this scenario, and (b) so our results are less opaque and the first-born BH can only attain spin through accretion following BH formation, and the second-born BH can only attain spin through BH–WR tidal spin-up. Although subdominant, this channel contributes up to ∼10% of our underlying population of merging BBHs in certain models, similar to the underlying populations of Neijssel et al. (2019). It does, however, contribute to a much smaller fraction of the detectable population (see, e.g., Figure 12 of Neijssel et al. 2019). We anticipate that proper treatment of this spin-up process would lead to spinning first-born and second-born BHs and large, positive effective spins. However, in our models, this channel favors the formation of near-equal-mass binaries, with ∼95% of systems having mass ratios of q > 0.75. Therefore, although the inclusion and proper treatment of this channel would affect the formation of near-equal-mass systems with high effective spins, we find it to be inefficient at generating systems with spinning primaries and significantly asymmetric masses.

Another evolutionary scenario not explored in detail here that can induce the spins in BH progenitors is through MT that occurs during the main sequence. In this Case A MT scenario, the primary star is in a tight binary with an orbital period at ZAMS of ∼0.5−1.2 day and is tidally locked when it overfills its Roche lobe on the main sequence. This strips the primary star of its envelope and causes it to become a fast-spinning WR star without any significant expansion. Although this channel has been shown to lead to the formation of high-mass X-ray binaries with appreciable BH spin, it will typically not lead to a merging BBH (Belczynski et al. 2012; Qin et al. 2019; Neijssel et al. 2021).

Throughout this work, we have assumed that AM transport in massive stars is highly efficient, leading to BHs formed from progenitors that do not undergo substantial tidal interaction to have near-zero spins. Although the numerical simulations of AM transport in high-mass stars through Taylor–Spruit dynamo and the majority of spins in BBH systems being small hint at efficient AM transport, the astroseismic measurements that encode this process are typically for low-mass stars. Furthermore, efficient AM transport may still lead to BH spins that are low but nonnegligible (a ≃ 0.1, Belczynski et al. 2020), and processes such as supernova fallback may induce some spin on the resultant BH even if AM transport is highly efficient (Schrøder et al. 2018). As a simple approximation for less efficient AM transport, we set a floor on the minimum spin of BHs at birth:

where a_min is the minimum spin magnitude of quasi-isolated BHs at birth. With a_min = 0.2, we find only minor changes in our main results regarding MRR and spinning primaries. For example, with all second-born BHs satisfying the criteria a_2b ≥ 0.2, we find that the percentage of systems in our most optimistic MT efficiency models with q_b > 2 and a_2b ≥ 0.2 to raise from ∼0.6%–1.5% to ∼0.7%–1.9%. However, setting the minimum natal spin magnitude of BHs has a larger impact on the systems that do not proceed through a MRR, with upwards of 52% of systems satisfying q_b < 0.5 and a_1b ≥ 0.2.

Of course, isolated evolution is only one of many proposed formation scenarios for generating BBH mergers, and the full population of BBH mergers may be the result of a combination of formation pathways (Bouffanais et al. 2021a; Wong et al. 2021; Zevin et al. 2021). Although we will not cover alternative potential formation pathways in detail here (see, e.g., Mandel & Farmer 2022 for a review), we note that the formation of systems with spinning primaries and asymmetric masses may be the result of hierarchical mergers in dense stellar environments. However, the formation of GW190412-like systems still proves to be enigmatic in the standard hierarchical merger paradigm, as its primary BH component is unlikely to be the product of a second-generation merger in an environment such as a classical globular cluster because its primary spin is lower than what one would expect for the merger product of two low-spinning BHs (Abbott et al. 2020). A number of other formation scenarios have been proposed such as third-generation mergers (Rodriguez et al. 2020), hierarchical systems (Hamers & Safarzadeh 2020), and formation in the disks of active galactic nuclei (McKernan et al. 2021; Tagawa et al. 2021). If systems such as GW190412 and GW190517_055101 result from subdominant formation scenarios, it may impact the astrophysical interpretation of the apparent anticorrelation between mass ratio and effective spin, as the exclusion of these events in such analyses (slightly) reduces the statistical significance of such a correlation (Callister et al. 2021), and different formation pathways do not necessarily need to share correlations in their mass and spin parameters.

The population of BBH mergers observed via GW emission has already led to many astrophysical lessons. In addition to the properties of systems as a whole, the distribution of mass and spin across the two components of the merging binary is enlightening when considering formation scenarios and physical processes. As always, the hundreds of anticipated BBH detections in the upcoming observing runs of the LVK interferometer network will help to determine whether such asymmetric-mass, highly spinning systems are commonplace. Given the quantity and flexibility of the proposed formation scenarios for compact binary mergers, we stress the importance of determining the regions of parameter space that cannot be easily populated by astrophysical formation models in addition to regions that can. This will be equally as important for constraining the merger rates and uncertain physical processes that embed the myriad of potential formation pathways.

We thank Mario Spera for insightful discussion and a detailed read of this manuscript, Luke Kelley for providing processed Illustris-TNG data used in our population resampling, and Katie Breivik for useful conversations as well as pertinent developments to COSMIC. We also thank the anonymous referee whose comments improved this manuscript. Support for this work and M.Z. was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51474.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. S.S.B. is supported by a Swiss National Science Foundation professorship grant (project number PP00P2_176868; PI Fragos).

Software: Astropy (Robitaille et al. 2013; Price-Whelan et al. 2018); iPython (Pérez & Granger 2007); Matplotlib (Hunter 2007); NumPy (Oliphant 2006; Van Der Walt et al. 2011); Pandas (McKinney 2010); SciPy (Virtanen et al. 2020).

The observation of compact binary coalescences via GWs is prone to selection effects, most notably that more massive systems are more luminous and thus can be seen out to higher luminosity distances. Current ground-based detectors are most sensitive at frequencies of ∼100 Hz, and systems with redshifted total masses ≳500 M_⊙ become unobservable since they merge at frequencies below the seismic floor (e.g., Mehta et al. 2022). Mass ratios and spins also impact the observability of compact binary mergers, with the detectors being most sensitive to systems with near-equal masses and high, aligned spins. Selection effects must be accounted for in our population models to yield a fair comparison against the observed BBH population.

We use a semianalytic treatment to incorporate selection effects in our population models, which relies on precomputed detection probabilities over a grid of chirp masses , mass ratios q = m₂/m₁ with m₂ ≤ m₁, and redshifts z where m₁ and m₂ are the primary and secondary masses. We assume a three-detector network consisting of LIGO–Hanford, LIGO–Livingston, and Virgo operating at either midhighlatelow or design sensitivity (Abbott et al. 2018). Extrinsic parameters are sampled over, and systems with network signal-to-noise ratio of ρ_thresh = 10 are considered detectable (Nitz et al. 2020; Abbott et al. 2021b), such that the detection probability is given

where ψ_j are the sampled extrinsic parameters, i defines the detector, and is a Heaviside step function. Spin effects are ignored, as they have a minor effect on detection probabilities (Ng et al. 2018). The relative weight of each system in the population thus becomes

where is the differential comoving volume at the merger redshift and is the time dilation between the merger and detectors. Detection probabilities for each system are determined using a k-nearest neighbors regressor trained on our detection probability grid.

Figure 8 shows the sensitive comoving spacetime volume, VT, for which a given search is sensitive to a BBH system with a particular primary mass and mass ratio, defined as

where T is the observing time, which is fixed to 1 yr, and is the detection probability for a system merging at redshift z with primary mass m₁ and mass ratio q (see, e.g., Fishbach & Holz 2017). Our grid of precomputed detection probabilities is publicly available on Zenodo (Zevin 2021), and the associated codebase for applying the selection effects is available on Github. ¹⁰

Zoom In Zoom Out Reset image size

Throughout this work, we have assumed that the critical mass ratios that govern the onset of unstable MT follow the prescription in Neijssel et al. (2019). Here, we briefly discuss how key results change with an alternative choice of critical mass ratios by adopting the values from Belczynski et al. (2008). Figure 9 shows the birth mass ratio versus effective spin distribution across our five simulated accretion efficiency (f_acc) values, for a single model variation of CE efficiency (α_CE = 1) and super-Eddington accretion factor (γ_Edd = 1). Overall, the mass ratio and effective spin distributions are similar between the two critical mass ratio models. The critical mass ratio parameterization of Belczynski et al. (2008) leads to slightly less systems proceeding through MRR and tidal spin-up. Compared to the critical mass ratios of Neijssel et al. (2019), at f_acc = 0.5, the Belczynski et al. (2008) model decreases the percentage of systems with q_b > 1 from 47% to 22%. The fraction of systems that proceed through MRR and have second-born spins of a_2b > 0.2 decreases by about a factor of 2.

Zoom In Zoom Out Reset image size