© The Authors 2024

1. Introduction

After the detection of binary black hole (BH) and binary neutron star (NS) mergers (Abbott et al. 2016, 2017a) by the ground-based gravitational-wave (GW) observatories LIGO (LIGO Scientific Collaboration 2015) and Virgo (Acernese et al. 2015), NS – BH (NSBH) systems had remained elusive until the detection of the NSBH merger events GW200115 and GW200105 (Abbott et al. 2021b). For GW200105, with a high-spin prior for the NS, the BH and NS masses are estimated to be $M_{\mathrm{BH}}=8.9_{-1.5}^{+1.2}{M}_{\odot }$ and $M_{\mathrm{NS}}=1.9_{-0.2}^{+0.3}{M}_{\odot }$, respectively. The BH spin magnitude is constrained as ${\chi }_{\mathrm{BH}}=0.{08}_{-0.08}^{+0.22}$ and its effective inspiral spin parameter ${\chi }_{\mathrm{eff}}=-0.{01}_{-0.15}^{+0.11}$. For GW200115, the binary mass components are estimated to be $M_{\mathrm{BH}}=5.7_{-2.1}^{+1.8}{M}_{\odot }$ and $M_{\mathrm{NS}}=1.5_{-0.3}^{+0.7}{M}_{\odot }$. The primary spin magnitude, ${\chi }_{\mathrm{BH}}=0.{33}_{-0.29}^{+0.48}$, is not well constrained, while the effective inspiral spin parameter is found to be $-0.{19}_{-0.35}^{+0.23}$. The latter indicates a possible misalignment between the BH spin and the orbital angular momentum. The misalignment angle is estimated to be $2.{30}_{-1.18}^{+0.59}\mathrm{rad}$ (Fragione et al. 2021). Mandel & Smith (2021) then reanalyzed the GW200115 signal with their astrophysically motivated priors and suggested a non-spinning BH for GW200115.

With the release of the third Gravitational-Wave Transient Catalog (GWTC-3), in addition to GW200105 and GW200115, NSBH merger candidates include GW190426 (Abbott et al. 2021a), GW190917 (Abbott et al. 2021c), and GW191219 (Abbott et al. 2023). GW190426 and GW190917 are marginal candidates because their probability of astrophysical origin, p_astro, is below 0.5, and GW191219 exhibits large model-dependent uncertainties in p_astro. In the GWTC-3 analysis, GW200105 was also listed as a marginal candidate because its p_astro is below 0.5. Nevertheless, it is still an interesting NSBH candidate with a clear outlier from the noise background (Abbott et al. 2023). The properties of these NSBH candidates are listed in Table 1. Apart from these events, GW190814 (Abbott et al. 2020; Huang et al. 2020; Zhou et al. 2021) and GW200210 (Abbott et al. 2023) are suspected to be NSBH mergers. Both systems possess a compact object that is within the mass range of 2 − 3 M_⊙ and close to the boundary between the heaviest NS and the lightest BH. In the near future, more NSBHs are expected to be detected in the fourth (O4) and fifth (O5) observing runs of ground-based GW detectors at improved sensitivities (Abbott et al. 2018), as well as, in the longer term, the next generation of ground-based observatories, such as the Einstein Telescope (Punturo et al. 2010; Hild et al. 2011) and the Cosmic Explorer (Abbott et al. 2017b; Reitze et al. 2019) and future space-based observatories, such as the Laser interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Thorpe et al. 2019), TianQin (Luo et al. 2016; Mei et al. 2021), and Taiji (Ruan et al. 2020).

Neutron star – black hole inspiral can present different GW signals, depending on whether the NS is tidally disrupted and whether the mass shedding happens before the NS crosses the innermost-stable circular orbit (ISCO; Kyutoku et al. 2010, 2011; Foucart et al. 2013; Tiwari et al. 2021). If the NS is not plunging as a whole into the BH but is tidally disrupted outside ISCO, the NSBH mergers may have potential electromagnetic counterparts (EMCs) such as kilonovae (Li & Paczyński 1998; Kawaguchi et al. 2020; Zhu et al. 2021; Darbha et al. 2021), short gamma-ray burst emission (Nakar 2007; Pannarale & Ohme 2014; Gompertz et al. 2020), long gamma-ray burst emission (Gottlieb et al. 2023), and radio emission (Piran et al. 2013; Hotokezaka et al. 2016). Generally, stiff NS equations of state corresponding to large NS radii, low-mass BHs, and high BH spins, which affect the position of the ISCO, increase the possibility of generating EMCs (Pannarale et al. 2011; Foucart et al. 2018; Bhattacharya et al. 2019; Román-Garza et al. 2021). The recent study of Fragione (2021) showed that only if the BHs are spinning fast and NSs have stiff equations of state would one expect a significant fraction of NSBH mergers to be associated with EMCs. However, in the classic isolated binary formation channel, it is commonly thought that if the BH forms first, it is likely to have a small spin (Fragos & McClintock 2015; Qin et al. 2018; Fuller & Ma 2019; Bavera et al. 2023) because the angular momentum of the core can be efficiently transported to the outer layer, which would be removed by mass transfer or winds. Alternatively, NSBHs where the NS forms first, which allows for a tidally spun-up BH progenitor, become candidates of EMC sources (Román-Garza et al. 2021; Hu et al. 2022). Despite the possibility, this corresponds to a very small portion of the whole NSBH population (Shao & Li 2021; Broekgaarden et al. 2021; Chattopadhyay et al. 2022).

Regarding the formation mechanisms of merging NSBH binaries, multiple channels have been proposed. These channels include the evolution of isolated field binaries; dynamical interactions in global clusters (Clausen et al. 2013; Ye et al. 2020), nuclear star clusters (Petrovich & Antonini 2017; Arca Sedda 2020; Wang et al. 2021), and young star clusters (Ziosi et al. 2014; Fragione & Banerjee 2020; Rastello et al. 2020; Santoliquido et al. 2020; Arca Sedda 2021); evolution from triples (Fragione & Loeb 2019; Trani et al. 2022); and chemical homogeneous evolution (Marchant et al. 2017). The isolated binary evolution channel is widely explored through binary population synthesis (BPS) methods (Tutukov & Yungelson 1993; Fryer et al. 1999; Voss & Tauris 2003; Dominik et al. 2012, 2015; de Mink & Belczynski 2015; Ablimit & Maeda 2018; Giacobbo & Mapelli 2018; Kruckow et al. 2018; Neijssel et al. 2019; Belczynski et al. 2020; Tang et al. 2020; Zevin et al. 2020; Broekgaarden et al. 2021; Shao & Li 2021; Iorio et al. 2023). These methods estimate a merger rate density of NSBHs consistent with the observations, given appropriate physical assumptions. However, relying solely on the merger rate density is not sufficient for constraining model uncertainties. The merger rate density is subject to many uncertainties in binary and stellar evolution such as common envelope (CE) evolution, mass transfer processes, supernova (SN) mechanisms, and natal kicks, as well as the metallicity-specific star formation rate density and its redshift evolution (see, e.g., Broekgaarden et al. 2021, for a review). Varying model assumptions can make the rate predictions plausible, but this alone is not enough to obtain decisive constraints on the uncertain physical processes without detailed modeling and systematic joint investigations of all GW observables. As a result, the relatively well measured GW observables such as chirp masses and effective spins, which include the information of individual compact object masses, BH spins, and spin-orbit tilt angles, provide valuable constraints on the astrophysical models.

In this work, we present detailed binary evolution models that follow the entire evolution of binaries from zero-age main sequence (ZAMS) to the formation of merging NSBHs for the first time. We used the population synthesis framework POSYDON (Fragos et al. 2023) to carry out a BPS study of NSBH systems that will merge within a Hubble time. We focused on NSBH populations at solar metallicity and investigated their formation channels in detail. With POSYDON, we can evolve binaries using extensive simulation grids of stellar structure and binary evolution simulations computed with the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). The detailed simulations allow us to track the changes of the internal structure of the stars, taking into account both single-star evolution and binary interactions, as well as the angular momentum transport in the interior of the stars and between binary components. This enables us to make a self-consistent estimation of mass-transfer rate and thus an accurate appraisal of the mass transfer stability. Furthermore, the feedback of mass changes from stellar winds or mass transfer on the stellar structure is modeled self-consistently (see, e.g., Bavera et al. 2023). Such model advancement allows us to infer the end of the mass-transfer phase more accurately, resolving potentially partially stripped envelopes (Klencki et al. 2022). When the binary ensues into a CE phase, we are able to calculate the binding energy self-consistently from the internal structure of the actual donor star model. Finally, we can estimate the properties of compact objects such as masses and BH spins based on the internal stellar-structure profiles of the immediate progenitor stars.

This paper is organized as follows. In Sect. 2, we briefly introduce POSYDON and describe the most important input physics in the code. We present individual examples of binaries leading to the formation of NSBH, as well as the population properties of merging NSBH in Sect. 3, and we show the impact of a set of model uncertainties in CE evolution and core-collapse prescriptions in Sect. 4. Finally, we discuss our findings in Sect. 5 and summarize them in Sect. 6.

2. The binary population synthesis code – POSYDON

We used the newly developed BPS code POSYDON to study merging NSBHs. In contrast to most rapid BPS codes, POSYDON evolves binaries based on detailed stellar and binary models. For a detailed description of the code, we invite the reader to consult Fragos et al. (2023). In what follows, we only summarize some of its key aspects.

POSYDON contains five pre-calculated stellar and binary evolution grids including two single-star grids -a hydrogen (H)-rich star grid and helium (He)-rich star grid- as well as three binary grids in circular orbits with various initial orbital periods, primary masses, and mass ratios; these are the hydrogen ZAMS binary grid; an H-rich star at the onset of Roche-lobe overflow (RLOF) and a compact object grid; and an He-rich ZAMS and a compact object (CO-HeMS) grid. The single-star grids are used to treat the non-interacting binaries in eccentric orbits (see Sect. 8.1 in Fragos et al. 2023). The binary grids were first post-processed so that the data size of the original binary evolution calculations is reduced and additional quantities used in BPS are computed. We then applied classification and interpolation algorithms to the outputs of the three grids, which enabled us to effectively interpolate between the pre-computed binary evolution simulations and estimate the evolution of any arbitrary binary. These algorithms currently allow for the interpolation between the initial and the final state of a binary evolution calculation. As a simpler alternative, we also provide the functionality to evolve individual binaries using nearest-neighbor matching. While the former approach is the default in population synthesis calculations, nearest-neighbor matching can be used in situations where for a specific binary we need to demonstrate its detailed time evolution (see, e.g., Figs. 1 and 2).

Fig. 1. Evolutionary path, from ZAMS to the merger, of an NSBH binary from channel Ia. From top to bottom, we show the evolution of the orbital period Porb, eccentricity e, component masses M1, 2, star radii R1, 2, surface angular velocity over the critical value (ω/ωcrit)1, 2, and the spin of the primary χ1. The solid blue line represents the star that evolves to form a BH, while the red dashed-dotted line represents the star that becomes an NS. In the panel of stellar mass evolution, we indicate the onset of the stable mass transfer (SMT), the first core collapse (CC1), the common envelope (CE) phase, and the second core collapse (CC2) with the short vertical lines.

Fig. 2. Same as Fig. 1, but for an NSBH binary formed through channel IIa.

2.1. Stellar physics

We summarize the most important stellar physics model assumptions including stellar winds and mixing adopted in our stellar models. Regarding stellar winds, we used the Dutch scheme in MESA for stars initially more massive than 8 M_⊙. Specifically, we adopted the wind prescription of de Jager et al. (1988) for cool stars with effective temperatures below 10 000 K and Vink et al. (2001) for hot stars with effective temperatures above 11 000 K. Between these two temperatures, we linearly interpolated the wind-loss rates. If the surface H fraction of the hot star was below 0.4, we employed the Wolf–Rayet wind mass-loss rates of Nugis & Lamers (2000). For stars initially below 8 M_⊙, the Dutch scheme was used again when their effective temperature exceeded 12 000 K. On the other hand, for stars cooler than 8000 K, Reimers (1975) wind with a scaling factor η_R = 0.1 was adopted for those on the red giant branch and the Bloecker (1995) wind with scaling factor η_B = 0.1 was adopted for those in the thermally pulsating phase. Between 8000 K and 12 000 K, we linearly interpolated the wind-loss rates between the cool and hot winds.

To treat the convective mixing, we adopted the mixing-length theory (MLT, Böhm-Vitense 1958) with a mixing length parameter of α = 1.93. For the overshooting mixing, we used the exponential decay formalism (Herwig 2000; Paxton et al. 2011) with the free parameter f_ov = 0.016 for stars below 4 M_⊙ (Choi et al. 2016) and f_ov = 0.0415 for stars more massive than 8 M_⊙. The latter is converted from the result of Brott et al. (2011), that is, α_ov = 0.335 via the rough ratio ∼10 between the step overshoot formalism parameter α_ov and the exponential decay formalism parameter f_ov (Claret & Torres 2017). We note that we backed up 0.008 times the scale height from the edge of the convective zone when assigning f_ov. We made a smooth transition between the two values of f_ov for stars within 4 − 8 M_⊙. Additionally, we adopted the Spruit–Tayler dynamo (Spruit 2002) to account for the efficient coupling between the stellar core and envelope. We treated rotation mixing and angular momentum transport following the MIST project (Choi et al. 2016).

2.2. Tidal interactions and binary mass transfer

Tidal interactions in a binary system affect the orbital evolution and the rotational properties of the stars. In a close detached binary, tidal forces efficiently circularize the orbit and lead to the synchronization between the stars’ rotation rates and the orbit. We followed Hut (1981) to calculate the evolution of binary separation, eccentricity, and star rotation, considering the difference between convective and radiative envelopes (Hurley et al. 2002; Qin et al. 2018). For the initial state of the binary models, we assumed the spins of two ZAMS stars are synchronized with the orbits.

By virtue of detailed simulations, we were able to model the mass transfer in a self-consistent way, considering the feedback of mass loss or gain along with the angular momentum change for both donors and accretors. For a main-sequence star, we adopted the contact scheme in MESA to calculate the mass-transfer rate during RLOF. In contrast, for an evolved star, when its center H abundance is below 10⁻⁶, we used the Kolb scheme (Kolb & Ritter 1990) to allow for a better treatment of the extended envelope. Based on the nature of the accretor, we treated the mass-accretion process differently. In the case of accretion onto a compact object, we limited the accretion at the rate corresponding to the Eddington luminosity, L_Edd, which can be written as

$$\begin{array}{c}\hfill L_{\mathrm{Edd}}=\frac{4\pi G{M}_{\mathrm{acc}}c}{\kappa },\end{array}$$(1)

where G is the gravitational constant, c is the speed of light, M_acc is the mass of the accretor, and κ is the opacity. The excess mass transferred is assumed to be lost as isotropic wind, carrying with it the specific orbital angular momentum of the accretor. In the case of accretion onto a non-degenerate star, we considered a rotation-dependent accretion. The transferred material carrying angular momentum will spin up the outer layer of the accretor, considering both ballistic and Keplerian disk mass transfer (de Mink et al. 2013). The accretor’s rotation will reduce the accretion rate by enhancing its stellar wind (Langer 1998). Most importantly, accretion is ceased when the accretor reaches critical rotation, implicitly boosting its winds as much as needed in order to not exceed ω/ω_crit ∼ 1.

2.3. Common-envelope evolution

Common-envelope evolution is crucial for the formation of close double compact objects. If the binary experiences unstable mass transfer and the donor is not a main-sequence star or a stripped He star, the binary will enter a CE phase. Then, the companion transports orbital energy and angular momentum to the CE, causing dramatic orbital shrinkage. If the CE is successfully ejected during this process, the binary will avoid a merger and go into a detached phase with the core of the donor and the companion. In our detailed binary evolution simulations, we considered the following conditions as the criteria for the onset of unstable mass transfer. Firstly, whenever the mass-transfer rate exceeded 0.1 M_⊙ yr⁻¹, we assumed it becomes unstable as the mass transfer is only expected to be more rapid after reaching this limit. Secondly, we checked if the star had mass loss from the second Lagrangian point (L2). The mass loss from the L2 point is considered to take away a significant amount of angular momentum from the binary leading to a rapid orbital shrinkage and triggering CE evolution (Tylenda et al. 2011; Nandez et al. 2014). For the case where the donor is a post-main-sequence star, we checked L2 RLOF following the implementation of Misra et al. (2020), while for the case of two main-sequence stars, we calculated the L2 radius using the prescription from Marchant et al. (2016). Thirdly, for a compact object accretor, if its photon trapping radius (Begelman 1979) extended beyond the Roche-lobe radius, the binary would enter CE evolution. Finally, if both binary components filled their Roche lobe, that is the binary was in a contact phase, and at least one of the stars had evolved beyond its main sequence, we assumed the binary had entered a CE phase.

To treat CE evolution, we used the α − λ formalism (Webbink 1984; Livio & Soker 1988). The value of λ_CE depends on the mass and evolutionary stage of the stars, affected by the internal structure and internal energy of the star (de Kool 1990; Dewi & Tauris 2000). At the onset of the CE phase, we calculated the value of λ self-consistently using the stellar profile. To encompass the model uncertainty on the boundary of the core and the envelope, we implemented three hydrogen abundances for the definitions of the He core-envelope boundary, which are X_H = 0.01, 0.1 and 0.3. We took the value 0.1 as the default. The α_CE parameter determines how efficient the orbital energy can be converted to eject the envelope. The value of α_CE is still uncertain and not necessarily constant for stars with different properties. We took α_CE = 1.0 as the default value. With these two parameters, we can calculate the binding energy of the envelope and the post-CE orbit where the transferred orbital energy is enough to disperse the envelope. If in the post-CE orbit computed by the α − λ formalism the core of the donor and the companion were not filling their Roche lobe, we considered the CE phase to be successful.

2.4. Core-collapse physics and compact-object formation

We adopted the Fryer et al. (2012) delayed core-collapse prescription by default to calculate the remnant baryonic mass of core-collapse SN (CCSN). The pre-SN carbon-oxygen core mass M_CO, core determines whether the star explodes and the fall-back fraction f_fb of the ejected mass. In addition, we explored other core-collapse prescriptions including the Fryer et al. (2012) rapid prescription and the Patton & Sukhbold (2020) prescription. The latter determines the explodability of a star based on the M_CO, core and the carbon mass fraction of the stellar core at the point of carbon ignition. For electron-capture SN (ECSN), we considered stars with the final He core masses in the range of 1.4 − 2.5 M_⊙ (Podsiadlowski et al. 2004) as their progenitors.

To covert the remnant baryonic mass to gravitational mass for the compact objects, we followed the method of Zevin et al. (2020). Furthermore, we assumed the maximum mass for an NS to be 2.5 M_⊙ (see discussion in Abbott et al. 2020, and references therein). Finally, we estimate the birth BH spins based on the final stellar profiles in our grids following the method described in Appendix D of Bavera et al. (2021), while we assumed the spins of NSs are zero for simplicity.

A newly born compact object may receive an SN natal kick because of asymmetric mass ejection and neutrino mass loss that can tilt the orbit, alter the orbital separation and eccentricity, or even disrupt the binary. By default, we adopted a natal-kick velocity following a Maxwellian distribution with a velocity dispersion of σ_CCSN = 265 km s⁻¹ (Hobbs et al. 2005) for CCSN and σ_ECSN = 20 km s⁻¹ for ECSN (Giacobbo & Mapelli 2019). The BH kick velocities follow the CCSN distribution, rescaled by a factor of 1.4 M_⊙/ M_BH. We calculated the orbital tilt θ introduced by the natal kick using Eq. (5) from Kalogera (1996), which is defined as the angle between the angular momentum of the post-SN orbit and the spin vector of the compact object.

After the formation of the BH, the secondary star may initiate a mass-transfer phase. In that case, we assumed that the secondary star’s spin would always be aligned with the first post-SN orbit. The accretion process onto BHs is expected to affect both the magnitude and the direction of the BH spins (Fragos et al. 2010). In addition, the accretion disk may become a warped disk due to the Lensing–Thirring effect (Scheuer & Feiler 1996), which would accelerate the alignment between the BH spin and the orbit. However, in the case of binaries with stellar-mass BHs, the effect can be very weak (King et al. 2005). We assumed that the transferred material has the specific angular momentum at ISCO. The BH’s angular momentum J can be divided into two components, along (J_y) and perpendicular (J_x) to the angular momentum of the orbit. For simplicity, we assumed the BH gains all the mass at once and the transferred angular momentum contributes to changing J_y. After accretion, the total angular momentum of the BH is equal to $\sqrt{J_{\mathit{y},\mathrm{acc}}^2+J_x^2}$, the tilt angle from the first SN, θ₁, becomes arccos(J_y, acc/J), and the updated BH spin χ_BH is determined by J. The simplified assumption might lead to the underestimation of the change on BH spin and tilt because the specific angular momentum of the ISCO is constantly changing during the accretion phase. However, we find that this has only a marginal impact on the majority of the binaries, as discussed in the next section.

When the second SN happens, the orbit is tilted again. The final tilt angle of the BH can be calculated from the two tilt angles θ₁ and θ₂:

$$\begin{array}{c}\hfill cos\theta =cos{\theta }_{1}cos{\theta }_2\pm sin{\theta }_{1}sin{\theta }_{2}cos\alpha ,\end{array}$$(2)

where α ∈ [0, π] is the angle of the relative position of two SN in the orbit, which we chose randomly as we attributed random mean anomaly for the two SN kicks. When α = 0 or π, the three vectors – the BH spin, the first post-SN orbital plane, and the second post-SN orbital plane – all lie in the same plane. This occurs either when both SNs take place at the same relative position within the orbit or at opposite sides of the orbit. In such cases, the final tilt is the addition or subtraction of the two tilt angles, which is also determined by whether the tilts are in the same or opposite directions.

With the information of component spins and orbital tilt angles, we can calculate the effective inspiral spin defined as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{eff}}=\frac{M_{\mathrm{BH}}{\chi }_{\mathrm{BH}}cos{\theta }_{\mathrm{BH}}+M_{\mathrm{NS}}{\chi }_{\mathrm{NS}}cos{\theta }_{\mathrm{NS}}}{M_{\mathrm{BH}}+M_{\mathrm{NS}}},\end{array}$$(3)

where θ_BH and θ_NS are the misalignment angles between spin vectors of compact objects and the vector along the final orbital angular momentum. θ_BH is equal to the final tilt angle θ in Eq. (2). Because we assume χ_NS = 0, χ_eff is determined by the component masses, BH spin, and the final tilt angle.

2.5. Initial binary properties

For each of our populations synthesis models, we evolved 10⁷ binaries, starting from two H-rich ZAMS stars. We randomly sample the initial masses of the primary stars, M_1, i in the range of 7 − 120 M_⊙, following a distribution of the initial mass function from Kroupa (2001). The secondary masses, M_2, i are in the range of 0.35 − 120 M_⊙ following a uniform distribution between the minimum and M_1, i (Kobulnicky & Fryer 2007). The distribution of the initial orbital separation is logarithmically flat between 1 and 10⁵ R_⊙, and we assumed zero eccentricity for the primordial binaries. To display the population properties of merging NSBHs, we applied a burst star formation history, which assumes all stars form at the same time.

3. Results

3.1. Formation channels

In our simulation, we found BHs always form first in merging NSBH systems. We identified two main formation channels for the merging NSBH systems. One involves stable mass transfer before the formation of the first-born compact object and CE evolution between a BH and a nondegenerate star (≈70%; the percentage refers to the number of merging NSBH systems per unit of simulated stellar mass) and the other one involves only a stable mass transfer episode(s) during the evolution (≈30%). Each channel can be further divided into two sub-channels based on the details of the mass transfer episodes, as we discuss in the following sections. For each channel, we illustrate the evolutionary path of a typical binary for the dominant sub-channel from a POSYDON simulation as an example.

3.1.1. Stable mass transfer and common envelope evolution channel

We refer to the channel involving stable mass transfer between two nondegenerate stars and a CE phase after the BH formation as channel I. The first stable mass-transfer phase occurs when the primary star begins expanding during the H shell-burning stage or the He core-burning stage. After the mass-transfer phase, the primary star loses its H-rich envelope through stellar winds. The stripped star then undergoes CCSN, forming a BH. After the first core collapse, the binary remains detached until the secondary star enters post-main-sequence phase and fills the Roche lobe. The second mass-transfer phase is unstable, leading to a CE phase. If the CE evolution is successful, the secondary star will lose the H-rich envelope with the ejection of the CE, leaving behind an He star. About 94% of the binaries undergo a third mass-transfer episode -a case-BB stable mass-transfer phase- between the BH and the He star, followed by the secondary star exploding and forming an NS. We refer to this channel as channel Ia. Alternatively, the remaining binaries evolve without undergoing any further mass-transfer episode after the CE phase. This is regarded as channel Ib. An illustration of the evolutionary phases are indicated at the top of Fig. 1, in which we show the evolutionary path of a typical NSBH merger evolving through channel Ia.

The binary in Fig. 1 initially has a primary star of 33.50 M_⊙, a secondary star of 23.45 M_⊙, and an orbital period of 72.0 d. The primary loses mass through stellar winds down to 28.50 M_⊙ before the mass-transfer phase initiates. The stable mass transfer starts at 5.795 Myr after the birth of the binary system, lasting for only ∼15 000 yr. The secondary is spun up quickly due to the accretion, gaining about 0.4 M_⊙ in total. The primary experiences CCSN at 6.257 Myr forming a BH of 10.98 M_⊙ with a negligible spin of χ_BH = 0.005, and the binary goes into the detached phase. As the secondary evolves, the mass transfer from the secondary onto the BH starts, shrinking the orbit. Soon after the onset of the mass transfer, the secondary radius exceeds the L2 Roche lobe radius, leading to a CE phase at 8.073 Myr. When the CE is ejected, the secondary is a stripped He-rich star of 9.97 M_⊙, and it expands at the end of the He-burning phase, leading to a case-BB mass-transfer phase. Eventually, the secondary collapses to a 2.20 M_⊙ NS at 8.766 Myr, and the binary merges at 552.3 Myr due to energy and angular momentum loss from GW radiation.

3.1.2. Stable mass transfer channel

We refer to the channel that involves only stable mass transfer episode(s) as channel II. About 96.4% of them undergo two stable mass-transfer phases. The binary first undergoes a stable case-A or case-B mass-transfer phase, resulting in the donors being partially stripped at the end of the main sequence. Consequently, the primary stars do not enter a giant phase as in channel I, but gradually shrink after evolving off the main sequence. After the formation of the BH, the binary experiences stable case-A or case-B mass transfer and some of them may experience both case-A and case-B mass transfer at the second mass-transfer phase. Finally, the secondary star explodes forming an NS. We refer to this channel as channel IIa. In cases where the primary stars are more massive than ≈60 M_⊙, the binaries may avoid the first mass transfer. This is because at solar metallicity, our POSYDON stellar models predict that the stars above ≈60 M_⊙ do not expand to the red supergiant state with radii of ∼1000 R_⊙ (Bavera et al. 2023). The mass loss through stellar winds from the stars keeps widening the orbits until the first core collapse. Subsequently, the binaries undergo case-A or case-B stable mass transfer from the secondary onto the BH and the secondary stars become NSs in the end. This sub-channel accounts for about 4.6% of channel II, and we refer to it as channel IIb. Without an efficient orbital shrinkage mechanism, at the time of NS formation, binaries from channel II have relatively wide orbits compared with those of channel I. Consequently, only the NSBHs formed with eccentric orbits (caused by SN natal kicks) result in a delay time shorter than a Hubble time.

Figure 2 shows a a typical binary evolving through the channel IIa, which initially has a 23.66 M_⊙ primary, an 11.83 M_⊙ secondary, and an orbital period of 13.90 d. The first stable mass-transfer phase starts at 7.867 Myr and ends at 7.882 Myr when the primary is still a main-sequence star. The secondary star reaches near-critical rotation during the mass-transfer phase, with the rate as high as ≈10^−2.5 M_⊙ yr⁻¹. After a short time of being detached from the secondary, the primary evolves off the main sequence, and the mass transfer starts again at 8.003 Myr and ends at 8.012 Myr with a peak of ≈10^−3.3 M_⊙ yr⁻¹. During the mass-transfer phase, the secondary star accretes about 0.4 M_⊙ in total. In the meantime, the primary loses its envelope through mass transfer and winds. At 8.600 Myr, the primary explodes to form a BH of 4.40 M_⊙ with a dimensionless spin parameter of χ_BH = 0.05. Then, a stable mass-transfer phase ensues from the secondary onto the BH. The secondary later explodes to form an NS of 1.28 M_⊙ at 21.73 Myr. Finally, the binary merges at 6116 Myr.

Our results regarding formation channels demonstrate disparities with previous population studies. In the investigations by Dominik et al. (2012) and Broekgaarden et al. (2021), which utilized rapid BPS codes based on the SSE fitting formulae (Hurley et al. 2000), the NSBH formation is heavily dominated by channel I (≳90%) at solar metallicity. Iorio et al. (2023), employing updated single stellar tracks in their simulations, predicted a reduced fraction of about 50% for channel I at solar metallicity. However, their fraction of channel II becomes marginal at high metallicities. These discrepancies stem from differences in stellar models implemented in BPS codes and the approach to treating mass transfer stability.

3.2. Population properties of merging NSBHs

In this section, we present the population properties of merging NSBHs and their progenitors, including component masses, orbital properties, BH spins, final tilt angles, and effective spins. We selected all the NSBH systems that merge within a Hubble time and distinguish them by the formation channels¹.

3.2.1. Progenitors of merging NSBHs

In the left panel of Fig. 3, we show the initial orbital periods P_orb, i versus the initial primary masses M_1, i of merging NSBHs. The majority of the binaries from channel I have P_orb, i above ≈50 d, and most of the binaries evolving through channel II have closer orbits with P_orb, i in the range from a few days to a few tens of days. Binaries with initial primary masses high enough to produce a BH and initial periods ≲50 d are more likely to result in post-BH-formation orbits that subsequently lead to stable case-A or case-B mass transfer, from the secondary star onto the BH this time. Even if the mass ratio is such that the mass transfer is unstable, the CE evolution will result in a merger, as the donor star is still very compact. For binaries with initial periods ≳50 d, even when they result in post-BH formation binaries that undergo stable mass transfer, the final orbital periods at the formation of NSBH that are too wide to lead to a merger within a Hubble time. In contrast, the progenitor binaries from channel I, must have initial periods ≳50 d, such that the post-BH formation orbital periods are wide enough to not only lead to unstable mass transfer, but they also have the secondary star sufficiently evolved at the onset of the unstable mass transfer so that CE evolution can lead to the successful ejection of the envelope. The interplay of all these processes determines the relatively sharp boundary in initial orbital periods that separate channel I and channel II.

Fig. 3. Initial properties of progenitors of merging NSBHs. The blue, green, orange, and red dots represent binaries from channels Ia, Ib, IIa, and IIb, respectively. The left panel shows the primary masses versus initial binary orbital periods. The gray line indicates a 50 day orbital period for reference. The right panel displays the component masses with the gray lines indicating different mass ratios q, defined as M2, i/M1, i.

In channel II, the binaries with the most massive primary stars have the shortest initial orbital periods of a few days. The strong winds of the massive stars will keep widening the orbit during the early evolution, resulting in orbital periods of ≈10 d in the end. These binaries experience no RLOF mass transfer before the BH formation because the massive primary stars do not expand to become supergiant stars.

The ZAMS masses of the progenitor stars of merging NSBHs are shown in the right panel of Fig. 3. The primary stars that evolve into BHs have initial masses M_1, i varying from ≈17 to ≈120 M_⊙, with the majority of them in the range of ≈20 to ≈40 M_⊙. The secondary stars that form NSs have initial masses M_i, 2 lying within the range of ≈8 to ≈30 M_⊙. Most of the binaries evolving through channel I exhibit initial mass ratios q = M_i, 2/M_i, 1 greater than 0.5, while binaries originating from channel II display a broader range of mass ratios. Notably, binaries from channel II have small mass ratios and can reach below 0.2. The right panel of Fig. 3 also shows the scarcity of massive (≳40 M_⊙) primary stars in channel I. Given the adopted stellar physics parameters in POSYDON, a 40 M_⊙ can lead to the formation of BH as massive as 18 M_⊙. Such massive BHs, when paired with secondary stars that are potential NS progenitors (≲28 M_⊙), they experience stable mass transfer independently of the post-BH formation orbital periods and thus cannot be progenitors of merging NSBHs from channel I.

3.2.2. Compact objects in merging NSBHs

The left panel of Fig. 4 shows the orbital periods and eccentricities of the merging NSBHs at the time of double compact object formation. Channel I is capable of forming merging NSBH systems in close obits ≲0.1 d due to CE evolution, whereas channel II can only form them in orbits ≳1 d. A large number of the NSBH systems forming without undergoing a CE phase have too long delay times to merge within a Hubble time, except for those that obtain a high eccentricity from natal kick at the second SN. The merging NSBHs from channel II have eccentricities ≳0.5.

Fig. 4. Properties of merging NSBHs during formation of double compact object system. The left panel shows the orbital periods and eccentricities. The right panel shows the BH and NS masses and the 90% confidence intervals of GW200115, GW200105, GW190917, and GW190426.

In Fig. 4 on the right, we show the BH and NS masses of merging NSBHs as well as the 90% confidence intervals of NSBH merger candidates GW200115, GW200105, GW190917, and GW190426. GW191219 is not included because its BH mass is much larger than the maximal values in our population. The median of the component masses of the four candidates are within the mass range of merging NSBHs from our model. However, we cannot rule out the possibility that they were formed in a lower metallicity environment. Overall, the NS masses cover a mass range from ≈1.26 to ≈2.4 M_⊙, and the BH masses are from ≈2.5 to ≈20 M_⊙. Channel I forms BHs less massive than ≈12 M_⊙, while channel II is able to produce more massive BHs, but only a small fraction of them are above 15 M_⊙. The Fryer et al. (2012) delayed mechanism we adopted to calculate the compact object mass from CCSN does not produce a BH mass gap from 2.5 − 5 M_⊙. The BHs in the mass gap account for 16.2% of the total NSBH population. In the Fryer et al. (2012) delayed prescription, the discontinuity of calculating proto-compact object mass at M_CO, core of 3.5 M_⊙ introduces an NS mass gap of ≈0.1 M_⊙ around 1.7 M_⊙. However, the existence of this gap is not observed in practice. In our simulation, the interpolation applied to obtain the NS masses smooths out the discontinuity. In the Fryer et al. (2012) delayed prescription, if the final M_CO, core of the compact object progenitor is less massive than 2.5 M_⊙, a constant NS mass of ≈1.27 M_⊙ is predicted. Accordingly, we can see the concentration at that mass especially for channel II. In channel I, the NS progenitors are more likely to have lower M_He, core in the range of 1.4 − 2.5 M_⊙, leading to ECSN and the formation of NSs at ≈1.26 M_⊙.

We calculated the BH spins using the stellar profiles of the BH progenitors at the onset of the SN event (see Sect. 8.3.4 in Fragos et al. 2023). The BH spins versus initial orbital periods is shown in Fig. 5, where the left panel shows the BH spins at birth and the right panel shows the final BH spins after accretion. From the left panel, we can see a significant correlation between the BH spins and the initial orbital periods where shorter orbital periods correspond to higher BH spins. The prevalence of short initial orbital periods among the binaries in channel II results in a substantial number of BHs with non-negligible spins. In our models, the BHs always form from the primary stars. As we assume efficient angular momentum transport inside stars in our models, the primary stars are expected to lose most of their angular momentum through mass transfer and winds, forming BHs with negligible spins. However, the orbits of the initially close binaries can remain tight before the first core collapse, and tides could help spin up the progenitor stars. In this case, although the stars lose the envelope with nearly all the stellar angular momentum, the stripped primary stars preserve a portion of angular momentum contributing to the BH spins. Nevertheless, neither channel produces BHs with spins larger than ≈0.2 at BH birth, consistent with (Qin et al. 2018; Fuller & Ma 2019; Belczynski et al. 2020).

Fig. 5. BH spins versus initial orbital periods of merging NSBHs at BH birth (left) and after NS formation (right).

From the right panel, we can see that BHs can be moderately spun up by accretion. For channel I, the BHs gain spins less than ≈0.3 through case-BB stable mass transfer. For channel II, the BH spins increase mainly due to case-A stable mass transfer. Most of the BH spins still stay below 0.2, 8.4% of them are within the 0.2 − 0.4 range and only 0.7% of them exceed 0.4. In channel II, the inclination angle between the BH spin and the orbit is typically small because the pre-SN orbit is tightly bound and hard to tilt significantly with natal kicks. In channel I, some binaries possess large tilt angles or are even anti-aligned after the BH formation. However, the BH spins are close to zero. As a result, the impact of the simplified assumption on the spin evolution due to accretion is negligible for most of the binaries.

3.2.3. Spin-orbit tilt and effective inspiral spin

In Fig. 6, we show the distribution of the tilt angles between the BH spin and the orbital angular momentum, and the effective inspiral spins for the two channels. The numbers shown in the figure are normalized with the total underlying simulated mass of the population in units of ${10}^{-7}{M}_{\odot }^{-1}$. A large fraction of binaries from channel I have small tilt angles, and only about 2% of them end up with anti-alignment between the BH spin vector and orbital angular momentum vector. In contrast, for channel II, about 40% of the binaries have final tilt angles > π/2 rad.

Fig. 6. Distribution of final tilt angles and effective inspiral spins for channel I (blue) and II (orange). The colored circles and corresponding error bars in the right panel show the median and 90% credible intervals of effective inspiral spins for the NSBH merger event candidates. For some of the observed systems, the error bars extend outside the range of the figure (see also Table 1).

The difference for the two channels is attributed to the different orbital velocity at the time of the second SN. In channel II, binaries experience a natal kick with a velocity comparable to the orbital velocity to provide high eccentricity to the orbit. Hence, the kick can tilt the orbit easily. In channel I, most binaries are in tightly bound orbits with high orbital velocities; thus, natal kicks in excess of ≳1000 km s⁻¹ would be required to tilt the orbit significantly. With the final tilt angles and BH spins, we calculated the effective spins of NSBH systems. The right panel of Fig. 6 shows the distribution of the effective inspiral spins for the two channels as well as the median and cut-off 90% credible intervals of effective inspiral spins for the NSBH merger event candidates. Most of the binaries from channel I possess BHs with nearly zero spins and small tilt angles, so the distribution of the effective spins peaks sharply at zero. A small fraction of BHs that have spun up due to accretion, contribute to the presence of systems with positive χ_eff ≳ 0.02. Channel II produces a larger fraction of spinning BHs and anti-aligned systems compared with channel I, leading to a flatter distribution of effective spins in the range of ≈ − 0.2 to ≈0.2. The number of the binaries exhibiting non-negligible positive χ_eff, for example, χ_eff ≳ 0.05, is comparable between the two channels. In contrast, it is only channel II that can produce NSBH mergers with χ_eff ≲ −0.02. Binaries with χ_eff below −0.1 and above 0.1 account for only ≈2% for channel I and ≈14% for channel II.

3.2.4. Associated electromagnetic counterparts

Observations of EMCs associated with NSBH mergers can provide unique insights into the GW physics and put constraints on the NS equation of state. To date, no EMC has been associated with identified NSBH merger candidates. To predict the fractions of NSBH mergers associated with an EMC in our simulations, we follow the method presented in Foucart et al. (2018), which involves calculating the remnant mass outside the ISCO. We considered three scenarios for the NS radii of 11 km, 12 km, and 13 km. In the calculations, we also took into account the misaligned angles between the BH spins and the orbits using the component of the BH spin that is aligned with the orbital angular momentum.

The first line of Table 2 shows the fractions of NSBH mergers that are associated with an EMC under the different assumptions of NS radii. When assuming an NS radius of R_NS = 11 km, the fraction of NSBH systems associated with an EMC for channel I is 1.69% and only 0.62% for channel II. In total, the fraction is 1.35%. These NSBHs feature light NSs < 1.29 M_⊙ and light BHs < 3.54 M_⊙. Even when assuming an NS radius of R_NS = 12 km, the BH masses in these NSBHs remain below 4.77 M_⊙. This suggests that the NSBH mergers associated with EMCs involve a BH that falls within the mass gap of 2.5 − 5 M_⊙. The largest NS radius assumption of R_NS = 13 km predicts a larger fraction of 15.80%, 9.83%, and 13.92% for channel I, channel II, and the total, respectively. The maximal BH masses of the binaries with EMCs reach up to 5.87 M_⊙, and the NS masses extend to 1.59 M_⊙. In this scenario, approximately 88%, 75%, and 85% of BHs in NSBH mergers associated with EMCs are within the mass gap for channel I, channel II, and the total, respectively.

Previous studies (Román-Garza et al. 2021; Hu et al. 2022) showed that if the NS forms first in an NSBH system, the BH progenitor can be spun up efficiently through tidal interactions with the NS, which results in a BH with high birth spin. The latter increases the possibility of having EMC. However, in our simulation, NSBHs with NS forming first account for a negligible fraction (< 0.01%). This is different from the results in Román-Garza et al. (2021), where the fraction of NSBHs with NS formed first can be higher than 10%. We attribute the difference to the modeling of accretion onto nondegenerate stars. While in Román-Garza et al. (2021) the rapid BPS code COSMIC (Breivik et al. 2020) was used to simulate the binaries from two ZAMS stars until the formation of compact object–stripped He star binaries, in this study, we used the detailed simulations throughout the evolution of the binaries. The mass accretion model adopted for nondegenerate accretors in Breivik et al. (2020), which only limits the accretion when it exceeds ten times the thermal-timescale mass-transfer rate of the accretor, is significantly more efficient compared with the rotation-dependent model (de Mink et al. 2013) adopted in POSYDON’s detailed binary evolution calculations. In this model, the accreted material carries with it the specific angular momentum of a Keplerian accretion disk whose inner edge is at the surface of the accreting star. As a star is accreting mass, it is being spun up due to the high specific angular momentum of the accreted material, and it typically reaches critical rotation after accreting only a fraction of a solar mass. After the accreting star reaches critical rotation at its surface, we assume that further accretion of material ceases. Overall, this process results in very inefficient mass transfer (≲10%), except in case-A mass transfer, where tidal spin down from the companion star counteracts the spin up due to accretion and accretion efficiencies of ≈30% can be reached (see also Langer et al. 2020).

The amount of mass that can be accreted onto the secondary, during the initial mass-transfer phase between two non-degenerate stars, plays a crucial role in determining whether the formation of the NS can precede the formation of the BH (Sipior et al. 2004; Broekgaarden et al. 2021). The accretion model adopted in POSYDON’s detailed binary evolution calculations results in low accretion efficiency and thus does not lead to enough mass accreted by the secondary stars to form BHs, while the primary stripped star forms an NS first. Shao & Li (2021) also predicted a very low fraction of NSBHs where the NS formed first, when using a similar rotation-dependent model. An additional factor that may lead to an increased fraction of first-born NSs in NSBH populations modeled by rapid codes is that case-A mass transfer, which is most often involved in the formation of NSBH binaries with first-born NSs, is not accurately modeled in rapid BPS codes (especially those using Hurley’s fitting formulae; Hurley et al. 2000). The primary reason for this is that the formulae do not track the mass of the developing He core of the main-sequence star and the radius evolution of the star used to infer the end of the mass transfer does not take into account altered structure of the stripped donor star. Both these factors can lead to more extended case-A mass-transfer cases, where a larger fraction of the donor star is transferred onto the accretor (see Dorozsmai & Toonen 2022, for an extended discussion).

In our detailed binary evolution calculations, the accretion efficiency in case-A mass-transfer phases between two nondegenerate stars can reach up to ≈30%, which is in principle sufficient to result in a first-born NS and a second-born BH. However, the orbits of those binaries after the formation of the NS are still compact, and the subsequent mass transfer from the secondary star onto the NS, which is unstable due to the high mass ratio, leads to a merger during the CE evolution.

4. Model variations

Binary evolution modeling is subject to various uncertainties that can affect the BPS results significantly. Thus, it is always instructive to investigate the effects of model variations. Because POSYDON allows on-the-fly calculations at the stage of core collapse and CE evolution, we are able to apply different models for core collapse and CE evolution, showing their impact on the population properties of merging NSBHs. In contrast, changing physical assumptions such as the accretion efficiency, as discussed in the previous section (see also Sect. 5.1), would require us to rerun the entire library of pre-calculated binary evolution models, and thus it is outside the scope of this work.

4.1. Variations in common envelope prescriptions

In CE evolution, we considered the variations of the core-envelope boundary and α_CE parameters. We explored two more values of the H abundance, which are 0.01 and 0.3 (the default value is 0.1), to define the core-envelope boundary. Different core-envelope boundaries will directly affect the binding energy of the envelope. With lower H abundance for core-envelope definition, the boundary will be deeper inwards the star and the binding energy of the envelope that needed to be overcome will increase, leading to a smaller λ_CE. In addition, as the default α_CE = 1.0, we apply a smaller α_CE parameter of 0.5 as a model variation for comparison.

In Fig. 7, we show the distributions of NS masses, BH masses, and chirp masses for merging NSBHs under different parameters associated with CE evolution. The variations in CE evolution do not affect the stable mass transfer process and thus overall distributions have no differences for channel II. We can see that the ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model produces BHs peaking at around 7 M_⊙, while the other two variations have a larger fraction of BHs around 10 M_⊙. The ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.01}$ model leads to a smaller core of the donor star and larger binding energy of the envelope, which makes the envelope harder to expel. A smaller value of α_CE increases the energy needed to expel the donor’s envelope. Consequently, a reduced α_CE effectively decreases the overall number of binaries from channel I. For these two models, a heavy BH is favored because it reserves more orbital energy in the same orbit. The ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model results in wider post-CE orbits, which increases the probability of having merger delay times longer than the Hubble time. In addition, the ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model produces high-mass progenitor stars for the NSs, which are reflected by the relatively small number of light NSs. Perhaps most notably, for the two model variations with less efficient CE evolution (models α_CE = 0.5 and ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.01}$), it is channel II that becomes the dominant one, with the fraction of merging NSBH binaries coming from channel I becoming 26.2% and 47.6%, respectively. Notably, as the model variations in CE do not impact channel II, the variability observed in each bin across the four models in channel II indicates the statistical uncertainties in our population. To assess this, we generated four random subsets for each model, estimating the statistical uncertainties to be consistent with assuming Poisson statistics errors. More importantly, we find that the characteristic features of the distributions are robust against statistical uncertainties.

Fig. 7. Distributions of BH masses, NS masses, and chirp masses of merging NSBHs, scaled by the stellar mass of the parent population, for different parameters associated with CE evolution. The top and bottom correspond to channels I and II, respectively. The default model corresponds to αCE = 1.0 and adopts ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.1}$.

In Fig. 8, we show the distributions of the final orbit tilt angle, effective spin, and delay time for different parameters associated with CE evolution. The variations in CE evolution have limited impact on the distribution of tilt angles and effective spins, other than the overall change in the normalization of the distribution. The merger delay time distribution, on the other hand, is sensitive to the α_CE parameter and core-envelope definition. A smaller α_CE or a smaller λ_CE means more orbital energy is needed to eject the envelope, suggesting a closer post-CE orbit, thus a shorter delay time. Under these two models, most binaries have delay time below ≈10^8.5 yr. On the contrary, we can see that the ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model predicts more merging NSBHs with long delay times ≳10^9.5 yr because the core-envelope boundary moves outward, the envelope’s binding energy decreases and CE evolution leads to less dramatic orbital shrinkage.

Fig. 8. Orbital tilt, effective spin, and delay time distributions for variations in CE evolution in channel I (top) and channel II (bottom).

We notice that with the ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model, merging NSBHs can also form through double CE evolution. These binaries have initial mass ratios close to unity and orbital periods of several tens of days. They will enter a contact phase after both stars have evolved off the main sequence, which we assume leads to a double CE phase. When the CE is ejected successfully, two stripped He cores will find themselves to be left in an orbit of a few days. Afterward, the two stars will remain detached till the formation of NSBHs. However, the number of the merging NSBHs originating from this channel is extremely limited; these only account for 0.1% of the total. We thus do not discuss this channel any further.

4.2. Variations in core-collapse prescriptions

The properties of the compact object formed from a star naturally depend on adopted prescription for core-collapse process. These prescriptions map the properties of the structure of the star right before the core collapse to the properties of the formed compact object, with different prescriptions predicting not only different compact object properties, but even different compact object types, for the same progenitor star structure. In this study, beside the default Fryer et al. (2012) delayed prescription, we considered two varying core-collapse prescriptions: the Fryer et al. (2012) rapid prescription and the Patton & Sukhbold (2020) prescription.

Similarly to in Fig. 7, we show the mass distributions for the two variation models compared with the default model in Fig. 9. The Fryer et al. (2012) rapid prescription produces a mass gap of 2.5 − 5 M_⊙ between NSs and BHs. It generates BHs within a mass range of ∼6 − 12 M_⊙ in channel I and an abundance of BHs with masses of ∼6 − 7 M_⊙ in channel II. The Patton & Sukhbold (2020) prescription predicts a BH mass range of ∼8 − 15 M_⊙ in both channels. Few NSs from the Fryer et al. (2012) rapid model exceed 2 M_⊙. For the Patton & Sukhbold (2020) prescription, the NS masses accumulate within 1.1 − 1.5 M_⊙. As for the chirp mass, the Fryer et al. (2012) delayed prescription predicts a distribution between 1.5 − 4 M_⊙. The Fryer et al. (2012) rapid prescription provides a narrower range between 2.5 − 4 M_⊙ because of the deficiency of the BHs in the mass gap and heavy NSs above 2 M_⊙. The Patton & Sukhbold (2020) prescription predicts a chirp mass range of 2.5 − 3.5 M_⊙ and a peak at ≈3 M_⊙.

Fig. 9. Same as Fig. 7, but with different core-collapse mechanisms. The default model corresponds to Fryer-delayed prescription.

Apart from the core-collapse prescription, a potential natal kick imparted onto the newly formed compact object also plays an important role in the formation of double compact objects. We considered two different values for the dispersion of the assumed Maxwellian kick velocity distribution for CCSN, σ_CCSN = 150 km s⁻¹ and σ_CCSN = 450 km s⁻¹ (the default value is σ_CCSN = 265 km s⁻¹). Figure 10 shows the tilt, effective spin, and delay time for the three different dispersion values of the CCSN natal-kick velocity distribution. Natal kicks significantly affect both the orbital properties of the NSBH systems and the overall merger rate. Smaller kicks decrease the possibility of binary disruptions or excessive delay times caused by SNs, thus enhancing the formation rate of merging NSBHs. Furthermore, small kicks produce a large number of merging NSBHs with small tilt angles. In channel II, compared with the default model, the model of σ_CCSN = 150 km s⁻¹ predicts more systems with aligned BHs, but almost the same number of anti-aligned systems. Larger kicks generate a flatter distribution of tilt angles, but the small total number limits the number of anti-aligned systems. The peak below 10⁸ yr in the delay time distribution corresponds to highly eccentric binaries. In contrast, we can see that in channel I, there is no significant secondary peak below short delay times (10⁸ yr). The binaries in these channels tend to have short orbits and high orbital velocities right before the formation of the second compact object. Hence, is it is very unlikely to form NSBH binaries with extreme eccentricities, which are required to achieve merger delay times of the order of 10⁷ yr. Kicks in excess of ≳1000 km s⁻¹ would be necessary for this scenario, which is why we only see a mild excess of very short merger delay times for the highest dispersion value of our assumed kick distribution.

Fig. 10. Same as Fig. 8, but with different natal kick velocity distributions for CCSN. The default model corresponds to σCCSN = 265 km s−1.

4.3. Fraction of associated electromagnetic counterparts

Table 2 lists the predicted fractions of NSBH mergers that could potentially produce EMCs under different models and for three assumed NS radii. We estimate the fractions in the range of 0 − 18.64%. The α_CE = 0.5 model and the ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.01}$ model yield lower fractions compared with the default model because in channel I, these two models predict a limited population of low-mass BHs (≲6 M_⊙). The ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.01}$ model also gives a smaller fraction, which can be attributed to a reduced proportion of low-mass NSs (≲1.5 M_⊙) in the model. The model of σ_CCSN = 150 km s⁻¹ predicts the highest fraction of NSBH mergers that are accompanied by an EMC. The variations of core-collapse prescriptions exert the most substantial impact, as they yield different mass distributions for the compact objects. The Fryer et al. (2012) rapid model does not produce BHs within the mass gap of 2.5 − 5 M_⊙. In the absence of low-mass BHs, the EMC fraction experiences a significant decline. Only when assuming the largest NS radius does the fraction in channel II become not marginal (5.08%). The BH masses are below 7.07 M_⊙ and their spins are in the range of 0.02 − 0.13. The inclination toward channel II is attributed to the low-mass BHs within this channel possessing slightly higher spins. The Patton & Sukhbold (2020) prescription predicts BH masses above ≈8 M_⊙, leading to negligible EMC fractions.

5. Discussion

5.1. Uncertainties in stellar and binary physics

Apart from the model variations in how we treat the core collapse and the CE phases discussed in the last section, other not fully understood stellar and binary evolution physics can also affect the BPS results in a significant way. The accretion efficiency onto a nondegenerate star is such an example for the formation of merging NSBH binaries. In our pre-computed grids of detailed binary evolution models, we adopted a rotation-dependent accretion model (de Mink et al. 2013), which assumes the accretion would be suppressed if the accretor reaches critical rotation. While this is a physically self-consistent way to treat the accretion of matter and angular momentum onto a star, some authors (Paczynski 1991; Popham & Narayan 1991) suggest that a critically rotating star can still accrete mass because the viscous coupling between the star and the accretion disk can regulate the angular momentum transport without affecting the mass accretion. In this case, the mass transfer efficiency can be much higher than our model. However, it is not clear whether the accretion would still be limited by other factors, such as the inflation of the accreting star due to high mass-transfer rate. Detailed binary evolution models that connect the accretion efficiency onto a star with such accretion disk models are not yet available.

The treatment of accretion onto a BH may also affect the properties of the NSBH population, especially the spin of the BH, which in turn determines whether the merger of the NSBH binary will be accompanied by an EMC. Radiation, general-relativistic, and magneto-hydrodynamics simulation of a supercritical accretion disk around BHs (e.g., Sądowski & Narayan 2016) suggest that the accretion rate may significantly exceed the Eddington limit. This possibility has been explored in the context of BH high-mass X-ray binaries (Moreno Méndez et al. 2008; Qin et al. 2022) and coalescing binary BH populations (e.g., Bavera et al. 2021; Briel et al. 2023). However, the effect of super-Eddington accretion on the formation and the properties of merging NSBH, and especially the probability of being accompanied by EMCs, is less well studied. We leave the exploration of different assumption regarding the accretion physics onto both nondegenerate stars and compact objects for a future study.

5.2. An additional formation channel

We found that in some binaries with an initial binary mass ratio close to unity, after a stable mass transfer from the primary to the secondary, the secondary evolves off the main sequence before the primary explodes and starts a reverse mass-transfer phase². At the beginning of the reverse mass-transfer phase, the accretor has been almost or totally stripped. It accretes hydrogen from the secondary, quickly creating a thin layer of hydrogen at the surface. The accretion tends to make the star spin up and expand. However, the accretor will reach critical rotation and cease further accretion before it has a chance to expand significantly. The remainder of the mass-transfer phase will proceed effectively as fully nonconservative mass transfer.

The behavior we described above depends completely on our assumption of mass-transfer efficiency on a nondegenerate star (see discussion in Sect. 5.1). If we consider instead the possibility that the accretion continues onto critically rotating stars, we would expect that the stripped He stars, after accreting a small amount of hydrogen, will swell up substantially and fill the Roche lobe. This would result in a CE evolution quickly after the onset of reverse mass transfer. Under the assumption that the reverse mass transfer would become unstable and initiate CE evolution immediately, we evolve these binaries further with POSYDON. It turned out that this evolutionary path can be important for the formation of merging NSBHs. The interesting part is that after the CE phase, two He stars will be left in a close orbit, and the two stars can be effectively spun up by tidal interactions. The BHs formed through this channel can have high spins at birth. Therefore, a thorough study on the accretion physics is needed in order to obtain the full picture of binary evolution and the formation of double compact objects.

5.3. Local merger rate density

In this study, we focused on models at solar metallicity. Thus, we did not calculate the merger rate density that can be directly compared with the estimate from GW detection. However, multiple works have shown that in contrast to BBH, the NSBH formation rate has limited connection with metallicity (Neijssel et al. 2019; Román-Garza et al. 2021; Broekgaarden et al. 2022; Iorio et al. 2023). In order to obtain a general idea on the local merger rate density from our simulation, we made the estimation by only integrating the star formation rate around solar metallicity. To do that, we followed the procedure described in Bavera et al. (2023) integrating the star formation rate in a metallicity range of [0.5 Z_⊙, 2 Z_⊙]. In this way, we underestimated the contribution of NSBH mergers at high redshifts where star formation at low metallicity is significant. Under these assumptions, our default model predicts a local merger rate density of 60.11 Gpc⁻³ yr⁻¹, while the model variations we considered generate local merger rate densities in the range of 17.00 Gpc⁻³ yr⁻¹ to 155.64 Gpc⁻³ yr⁻¹.

5.4. Observation of NSBH populations with electromagnetic waves

Before merging, the NS may appear as a radio pulsar, making the system a pulsar-BH binary. Detection of pulsar-BH binaries can give us valuable insights on the formation and evolution of NSBH binaries, extending our objects to different populations of NSBH binaries. BPS studies have estimated up to thousands of pulsar-BH binaries in the Milky Way (Chattopadhyay et al. 2021; Shao & Li 2021; Pol et al. 2021). They could be observed by the radio telescopes such as the Square Kilometre Array (SKA; Kramer et al. 2004), MeerKAT (Booth et al. 2009), the Five-Hundred Metre Aperture Spherical Radio Telescope (FAST; Nan et al. 2011), and the Arecibo radio telescope (Cordes et al. 2006). Pulsar-BH binaries with orbital periods shorter than a day may be observed by both radio telescopes and GW detector LISA (Chattopadhyay et al. 2021). These joint observations could be further used to constrain model uncertainties in the formation of NSBH systems.

6. Summary

In this study, we investigated the formation channels and population properties of merging NSBH systems formed from isolated binary evolution at solar metallicity using the newly developed BPS code POSYDON. POSYDON incorporates detailed stellar and binary models calculated based on the binary-star evolution code MESA. The main findings of our work are:

• We identify two main formation channels of NSBH mergers. In channel I (≈70%), the binaries undergo a stable mass transfer prior to the formation of the BH, a CE phase after the formation of the BH, and a subsequent case-BB stable mass-transfer phase (channel Ia, ≈94%) or no case-BB mass transfer (channel Ib, ≈6%). In channel II (≈30%), the binaries go through only stable mass transfer. Most binaries within channel II experience two stable mass-transfer phases, one prior to and another after the BH formation (channel IIa, ≈95.4%). In addition, a small portion of the binaries, which contain massive primary stars above ≈60 M_⊙, are stripped because of strong stellar winds and avoid mass transfer before the BH formation (channel IIb, ≈4.6%).

• Independently of the formation channels and the explored model uncertainties, our models suggest that the BHs always form first in merging NSBH systems. The primary reason for this difference compared with other NSBH BPS studies in the literature is our adopted model of rotation-dependent accretion efficiency onto nondegenerate stars, which generally results in highly nonconservative mass transfer.

• Most binaries in channel I have initial orbital periods ≳50 d and mass ratios ≳0.5. Given the initially large orbital separation of these binaries, the majority of the BHs form with near-zero spins, assuming efficient angular momentum transport within the stars. After their formation, BHs can still be spun up due to mass accretion during the case-BB mass-transfer phase. However, since we assume that the accretion onto a compact object is Eddington limited, this effect is not significant.

• Compared with channel I, most of the binaries formed through channel II have initial orbital periods ≲50 d and a wider spread of initial mass ratios. Those binaries end up in relatively wide orbits before the second core collapse, as the second stable mass-transfer phase does not shrink the orbit as efficiently as the CE. Natal kicks during the second SN can significantly alter those orbits and in some cases result in highly eccentric NSBH binaries. It is the high eccentricity (> 0.5) of the orbits after the NS formation that causes the binary to merge within a Hubble time. Furthermore, binaries in channel II, due to tides in their close initial orbits and the further possible shrinking of the orbit during the first mass-transfer episode, have primary stars that can retain some angular momentum until the end of their lives. The angular momentum preserved in the core contributes to nonzero small spins (≲0.2) for the BHs. The BHs can be further spun up moderately because of the mass-transfer phase following the BH formation.

• Binaries from channel I mostly have small orbital tilt angles, while for channel II, binaries with large tilt angles ( > π/2 rad) account for a significant fraction. The reason for the latter is that the NSBH binaries that experienced a sufficiently high kick to become highly eccentric manage to merge within a Hubble time in channel II, and the kick can easily tilt the orbit in such wide systems. This feature also contributes to the distinctive distributions of effective spin for the two channels. For channel I, the negligible BH spins and the small tilt angles make the distribution of effective spin peak at zero, with an excess at positive values. In contrast, for channel II, the effective spins extend to ≈ − 0.2 and ≈0.2.

• With the information of compact object masses, BH spins, and orbital tilt angles, we estimate the fraction of merging NSBHs that we expect to be accompanied by an EMC. We find that when assuming a stiff equation of state for the NSs (R_NS = 13 km), the EMC fraction is about 13.92%, while for a soft equation of state (R_NS = 11 km) it is about 1.35%. The binaries that have potential EMC favor light BHs and NSs. Approximately 85% of BHs are within the mass gap of 2.5 − 5 M_⊙ when assuming R_NS = 13 km, and in the other two cases all BHs are expected to fall within the mass gap. Incorporating uncertainties in both the population modeling and the neutron-star equation of state, we estimate that between 0 and 18.6% of NSBH mergers could potentially exhibit an EMC.

• In addition to our default BPS model, we investigated a range of model uncertainties, including different prescriptions and model parameters for the CE and core-collapse phase. Different parameters for CE evolution affect channel I in terms of the absolute number of merging NSBHs and the compact object masses. The X_H = 0.01 core-envelope boundary definition results in a smaller core and thus higher binding energy of the envelope compared with the default 0.1 X_H, while a smaller α_CE means more orbital energy is needed to overcome the binding energy. These two variations predict fewer merging NSBHs and a preference for massive BHs ≈10 M_⊙. The ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.3}$ model leads to a larger post-CE orbital separation. The formation rate of merging NSBHs is lower because of excessive delay time for binaries going through CE evolution. Variations in the core-collapse prescriptions affect both channels concerning the formation efficiency of merging NSBHs, as well as their mass distributions and orbital properties. With our default Fryer et al. (2012) delayed model, most merging NSBHs have BHs less massive than ≈11 M_⊙, chirp masses in a range of 1.5 − 4 M_⊙, and NS masses within 1.26 − 2.4 M_⊙. The Fryer et al. (2012) rapid core-collapse prescription produces the mass gap between NSs and BHs and NSs less massive than ≈2 M_⊙. The Patton & Sukhbold (2020) prescription produces a narrow NS mass range of 1.1 − 1.5 M_⊙. Assuming a Maxwellian distribution of natal-kick velocities with a smaller dispersion (150 km s⁻¹ versus the 265 km s⁻¹ of our default model) predicts more merging NSBHs with a larger proportion of them exhibiting small tilt angles. Besides, a smaller kick velocity dispersion produces more merging NSBHs with long delay times, as it is less efficient in forming very eccentric NSBH binaries.

The utilization of detailed stellar and binary models in BPS studies enables enhanced modeling accuracy and self-consistent treatment of stellar evolution and binary interactions, providing us with more accurate and informative results. As we expect a growing number of GW events to be detected in the future, it becomes imperative to employ increasingly accurate BPS models; this is a crucial component in the unraveling of the origin and formation mechanisms of observed double compact objects and in refining our understanding of underlying stellar and binary physics.

Acknowledgments

We thank Christopher Berry for helpful comments. The POSYDON project is supported primarily by two sources: the Swiss National Science Foundation (PI Fragos, project numbers PP00P2_211006 and CRSII5_213497) and the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477). Z.X. acknowledges support from the Chinese Scholarship Council (CSC). S.S.B., T.F., M.K., and Z.X. were supported by the project number PP00P2_211006. S.S.B. was also supported by the project number CRSII5_213497. A.D., K.A.R., P.M.S. and M.S. were supported by the project number GBMF8477. K.K. acknowledges support from the Spanish State Research Agency, through the María de Maeztu Program for Centers and Units of Excellence in R&D, No. CEX2020-001058-M. E.Z. acknowledges funding support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 772086). The computations were performed at Northwestern University on the Trident computer cluster (funded by the GBMF8477 award) and at the University of Geneva on the Yggdrasil computer cluster. This research was partly supported by the computational resources and staff contributions provided5for the Quest high-performance computing facility at Northwestern University, jointly supported by the Office of the Provost, the Office for Research and Northwestern University Information Technology.

References

All Tables

All Figures

Fig. 1. Evolutionary path, from ZAMS to the merger, of an NSBH binary from channel Ia. From top to bottom, we show the evolution of the orbital period Porb, eccentricity e, component masses M1, 2, star radii R1, 2, surface angular velocity over the critical value (ω/ωcrit)1, 2, and the spin of the primary χ1. The solid blue line represents the star that evolves to form a BH, while the red dashed-dotted line represents the star that becomes an NS. In the panel of stellar mass evolution, we indicate the onset of the stable mass transfer (SMT), the first core collapse (CC1), the common envelope (CE) phase, and the second core collapse (CC2) with the short vertical lines.

In the text

	Fig. 2. Same as Fig. 1, but for an NSBH binary formed through channel IIa.
In the text

Fig. 3. Initial properties of progenitors of merging NSBHs. The blue, green, orange, and red dots represent binaries from channels Ia, Ib, IIa, and IIb, respectively. The left panel shows the primary masses versus initial binary orbital periods. The gray line indicates a 50 day orbital period for reference. The right panel displays the component masses with the gray lines indicating different mass ratios q, defined as M2, i/M1, i.

In the text

	Fig. 4. Properties of merging NSBHs during formation of double compact object system. The left panel shows the orbital periods and eccentricities. The right panel shows the BH and NS masses and the 90% confidence intervals of GW200115, GW200105, GW190917, and GW190426.
In the text

	Fig. 5. BH spins versus initial orbital periods of merging NSBHs at BH birth (left) and after NS formation (right).
In the text

	Fig. 6. Distribution of final tilt angles and effective inspiral spins for channel I (blue) and II (orange). The colored circles and corresponding error bars in the right panel show the median and 90% credible intervals of effective inspiral spins for the NSBH merger event candidates. For some of the observed systems, the error bars extend outside the range of the figure (see also Table 1).
In the text

	Fig. 7. Distributions of BH masses, NS masses, and chirp masses of merging NSBHs, scaled by the stellar mass of the parent population, for different parameters associated with CE evolution. The top and bottom correspond to channels I and II, respectively. The default model corresponds to αCE = 1.0 and adopts ${\lambda }_{\text{CE}}^{X_{\text{H}}=0.1}$.
In the text

	Fig. 8. Orbital tilt, effective spin, and delay time distributions for variations in CE evolution in channel I (top) and channel II (bottom).
In the text

	Fig. 9. Same as Fig. 7, but with different core-collapse mechanisms. The default model corresponds to Fryer-delayed prescription.
In the text

	Fig. 10. Same as Fig. 8, but with different natal kick velocity distributions for CCSN. The default model corresponds to σCCSN = 265 km s−1.
In the text