Population Synthesis of Black Hole Binaries with Compact Star Companions

Stars are believed to end their lives as compact stars (CSs) that include white dwarfs (WDs), neutron stars (NSs) and black holes (BHs). Observations indicate that massive stars are predominately born in binary and multiple systems (Sana et al., 2012; Kobulnicky et al., 2014; Moe & Di Stefano, 2017). It is naturally expected that CS pairs can be formed through isolated binary evolution when both components have evolved to be CSs. Very close CS pairs emit gravitational wave (GW) signals that are likely to be identified by future space-based detectors e.g. LISA (Amaro-Seoane et al., 2017), TianQin (Luo et al., 2016) and Taiji (Ruan et al., 2020), and/or by ground-based interferometers e.g. LIGO (Aasi et al., 2015) and Virgo (Acernese et al., 2015). Common envelope (CE, see a review by Ivanova et al., 2013) evolution is thought to be a vital stage for the formation of close CS pairs, during which the binary systems rapidly shrink inside a shared envelope originating from the donor star and the orbital energy is dissipated to expel the CE (Paczynski, 1976; Webbink, 1984; Iben & Livio, 1993). Usually dynamically unstable mass transfer between binary components can result in the occurrence of a CE phase, which is critically dependent on the stellar properties of both the donor and the accretor, and the mass ratios of binary components (e.g., Soberman et al., 1997; Ge et al., 2010, 2020; Shao & Li, 2014; Pavlovskii & Ivanova, 2015; Pavlovskii et al., 2017). It is generally accepted that a binary system is likely to enter CE evolution if the mass ratio of binary components is extreme or the donor star has developed a deep convective envelope prior to the mass transfer.

Radio and X-ray observations of Galactic NSs suggest a maximum mass of $\sim 2M_{\odot}$ (Antoniadis et al., 2013; Alsing et al., 2018; Farr & Chatziioannou, 2020), which is consistent with the maximum stable NS mass inferred from observations of the NS–NS merger GW170817 (Margalit & Metzger, 2017; Ruiz et al., 2018; Shibata et al., 2019). Meanwhile, Galactic BHs in X-ray binaries are inferred to have a minimal mass of $\sim 5M_{\odot}$ (Bailyn et al., 1998; Özel et al., 2010; Farr et al., 2011). This led to the mass gap ($\sim 2-5M_{\odot}$) between NSs and BHs. However, recent evidence suggest that the mass gap is being populated from both electromagnetic (Thompson et al., 2019; Wyrzykowski & Mandel, 2020; Jayasinghe et al., 2021) and GW observations (GW190814, Abbott et al., 2020a). Whether there is the mass gap or not can shed light on supernova mechanisms for the formation of NSs and BHs. Although Gaia astrometry was already claimed to be able to discover invisible BHs with optical companions (e.g., Gould & Salim, 2002; Mashian & Loeb, 2017; Breivik et al., 2017; Yalinewich et al., 2018; Yamaguchi et al., 2018; Breivik et al., 2019; Shao & Li, 2019; Andrews et al., 2019; Wiktorowicz et al., 2020), merging BH–CS binaries that appear as GW sources are alternatively potential objects to test the mass gap and relevant supernova mechanisms.

There is no BH–WD binary identified so far, although they are proposed to be associated with many kinds of astrophysical objects, e.g. ultracompact X-ray binaries (UCXBs, van Haaften et al., 2012), subluminous type I supernovae (Metzger, 2012), tidal disruption events (Fragione et al., 2020), and long gamma-ray bursts (Dong et al., 2018). Some X-ray binaries in globular clusters are suggested to be BH–WD systems (Maccarone et al., 2007; Miller-Jones et al., 2015), and the formation of such systems in a dynamical environment has been explored by Ivanova et al. (2010). It is also noted that mass-transferring BH–WD binaries are potential GW sources that may be identified by LISA in the future (Sberna et al., 2021).

Since the discovery of the first GW signal from a BH–BH inspiral (Abbott et al., 2016), the LIGO/Virgo collaboration has reported the detection of tens of CS pair mergers (Abbott et al., 2019, 2020b). Among them the majority are confirmed to be BH–BH systems and several are suspected to be BH–NS systems, e.g. GW190814 (Huang et al., 2020; Zhou et al., 2020) and GW190425 (Han et al., 2020a; Kyutoku et al., 2020). As expected, searches of the electromagnetic counterparts for BH–BH mergers have yielded negative results (e.g., Copperwheat et al., 2016; Smartt et al., 2016; Racusin et al., 2017). It is predicted that LISA will detect a number of merging BH–BH binaries in the Milky Way, but that electromagnetic observations will be challenging (Lamberts et al., 2018).

The detection of BH–NS mergers is of great interest as they are expected to emit across a broad electromagnetic spectrum and have been suggested to produce radio flares (Nakar & Piran, 2011; Hotokezaka et al., 2016), kilonovae (Li & Paczyński, 1998; Zhu et al., 2020), short gamma-ray bursts (Paczynski, 1986; Gompertz et al., 2020), and so on. A recent quite comprehensive study of BH–NS mergers can be found in Broekgaarden et al. (2021). Before merger, Galactic BH–NS binaries are possibly observed as binary radio pulsars (Shao & Li, 2018). Close BH–NS systems in the Milky Way may emit GW signals with frequencies in the LISA band (Chattopadhyay et al., 2020).

With the growing population of GW sources, a large number of formation channels for merging CS pairs have been put forward, especially for the case of BH–BH systems. Until now, it is still impossible to distinguish between various channels based on GW observations alone. The popular channels involve the formation through isolated binary evolution (e.g., Tutukov & Yungelson, 1993; Lipunov et al., 1997; Voss & Tauris, 2003; Kinugawa et al., 2014; Belczynski et al., 2016; Eldridge & Stanway, 2016; Stevenson et al., 2017; van den Heuvel et al., 2017; Klencki et al., 2018; Mapelli & Giacobbo, 2018; Spera et al., 2019; Breivik et al., 2020; Zevin et al., 2020; Tanikawa et al., 2021; Bavera et al., 2021; Olejak et al., 2021), and dynamical interactions in globular clusters (Downing et al., 2010; Rodriguez et al., 2016; Askar et al., 2017; Perna et al., 2019; Kremer et al., 2020) or young stellar clusters (Ziosi et al., 2014; Di Carlo et al., 2019; Santoliquido et al., 2020; Rastello et al., 2020). Other channels include the formation from isolated multiple systems (Silsbee & Tremaine, 2017; Liu & Lai, 2018; Hoang et al., 2018; Fragione & Loeb, 2019), the chemically homogeneous evolution for rapidly rotating stars (Mandel & de Mink, 2016; de Mink & Mandel, 2016; Marchant et al., 2016), as well as the disk of active galactic nuclei (Antonini & Rasio, 2016; Stone et al., 2017; McKernan et al., 2018).

In this paper, we investigate the formation of merging BH–CS binaries via isolated binary evolution. Our previous works have focused on the Galactic populations of detached BH binaries with a normal-star companion (Shao & Li, 2019) and BH X-ray binaries (Shao & Li, 2020), by adopting the rapid (Fryer et al., 2012) and the failed supernova mechanisms (Sukhbold et al., 2016; Raithel et al., 2018) for NS/BH formation. Here we further incorporate the delayed (Fryer et al., 2012) and the stochastic (Mandel & Müller, 2020) recipes to deal with compact remnant masses and possible natal kicks, and then study the predicted properties of descendent BH binaries with a CS companion. For mass-transferring binaries with a BH accretor, Pavlovskii et al. (2017) demonstrated that the process of Roche lobe overflow (RLO) may be stable over a much wider parameter space than previously thought. Accordingly, Olejak et al. (2021) suggested that sufficiently strong constraints on mass transfer stability are necessary to draw fully reliable conclusions for the population of double CS mergers. However, considering that limited binary systems were included by Pavlovskii et al. (2017), in this work we will evolve a large grid of the initial parameters for BH binaries with a nondegenerate donor to deal with mass transfer stability and obtain thorough criteria (or parameter spaces) for the occurrence of CE evolution.

The structure of this paper is organized as follows. Section 2 describes the method of our binary population synthesis (BPS, see a review by Han et al., 2020b) and the input physics implemented in the BPS method, especially supernova mechanisms for NS and BH formation. In Section 3 we obtain the criteria of mass transfer stability for the BH binaries with detailed binary evolution calculation. Sections 4 and 5 show our BPS calculation results about the merging BH–CS binaries in the Milky Way and the local Universe, respectively. Finally we conclude in Section 6.

To simulate the formation and evolution of BH–CS binaries, we utilize the BSE code originally developed by Hurley et al. (2002) and significantly updated by Shao & Li (2014). Some further modifications on the code can be found in Shao & Li (2019) and Shao & Li (2020). We briefly summarize the most important points in the following. When dealing with the process of mass transfer in the primordial binaries, we use the rotation-dependent mode (Shao & Li, 2014) that assumes the accretion rate of the secondary stars to be dependent on their rotational velocities (see also Stancliffe & Eldridge, 2009). With this mode we are able to reproduce the distributions of known Galactic binaries including BH–Be star systems (Shao & Li, 2014, 2020) and Wolf-Rayet star–O-type star systems (Shao & Li, 2016). Because a large fraction of the transferred mass is lost from the rapidly rotating secondary star, the maximal mass ratio of the primary to the secondary stars for stable mass transfer can reach up to $\sim 6$, significantly larger than in the conservative mass transfer case (Shao & Li, 2014). During CE evolution, we employ the binding energy parameter $\lambda$ calculated by Xu & Li (2010) and set the CE ejection efficiency $\alpha_{\rm CE}$ to be unity. We follow Belczynski et al. (2010) to deal with the wind mass-loss rates for different types of stars, except that for helium stars we decrease the mass-loss rates of Hamann et al. (1995) by a factor of 2 (Kiel & Hurley, 2006).

For the formation of NSs and BHs, we take into account three supernova models to treat the remnant masses and natal kicks. These models involve (i) the rapid explosion mechanism (Fryer et al., 2012), (ii) the delayed explosion mechanism (Fryer et al., 2012) and (iii) the stochastic recipe developed by Mandel & Müller (2020). The rapid model predicts a dearth in the remnant masses between $\sim 2M_{\odot}$ and $\sim 5M_{\odot}$, while the other two models are able to produce CSs within the mass gap. For both the rapid and the delayed mechanisms, the remnant masses are determined by the CO core masses at the time of explosions, and subsequent accretion of the fallback material. We follow Fryer et al. (2012) to convert between the baryonic and gravitational masses for NSs (see also Timmes et al., 1996) and simply approximate the gravitational mass with $90\%$ of the baryonic mass for BHs. The maximum mass of NSs is set to be $2.5M_{\odot}$. We adopt the criterion of Fryer et al. (2012) to distinguish the NSs originating from core-collapse and electron-capture supernovae. NSs from core-collapse supernovae are assumed to be subject to a kick with a Maxwellian distribution with $\sigma=265\rm\,km\,s^{-1}$ (Hobbs et al., 2005), while NSs from electron-capture supernovae have a lower kick velocity with $\sigma=50\rm\,km\,s^{-1}$ (Dessart et al., 2006). For the natal kicks to newborn BHs, we use the NS kick velocities reduced by a factor of $(1-f_{\rm fb})$, where $f_{\rm fb}$ is the fraction of the fallback material. In the stochastic model, the outcome of supernova explosions is expected to be probabilistic rather than deterministic. The remnant masses and natal kicks are required to satisfy some specific probability distributions depending on the masses of the CO cores. Meanwhile, the hydrogen shell (if present) is assumed to be always ejected, and this will cap the BH masses at the helium core masses. This model allows a significant tail of low kicks for natal NSs, which can be consistent with the observation of large numbers of NSs in globular clusters (Mandel & Müller, 2020). In addition, a large fraction of BHs are expected to receive zero natal kicks. Following Mandel & Müller (2020), we take the maximum mass of NSs to be $2M_{\odot}$¹¹1The maximum mass of observed NSs is a bit higher, around $2.1M_{\odot}$, but this does not influence our final results..

Figure 1 shows the relations between the zero-age main sequence mass and the remnant mass for our adopted three supernova models. For the models involving the rapid and delayed explosion mechanisms, similar relations have been reported by many previous works (e.g., Belczynski et al., 2010; Fryer et al., 2012; Giacobbo & Mapelli, 2018; Shao & Li, 2019; Zevin et al., 2020). It is noted that the minimum mass of the BH’s progenitors can reach as low as $\sim 12M_{\odot}$ in the stochastic model. Since we follow Mandel & Müller (2020) to ignore the correction between the baryonic and gravitational masses, the stochastic model tends to produce slightly heavier BHs at the high mass end compared with those in the other two models. In our BPS calculations, we identify mass-gap BHs with masses between $2M_{\odot}$ and $5M_{\odot}$ for all models.

The evolutionary fate of the primordial binaries is determined by their initial parameters: the primary masses $M_{\rm 1}$, the secondary masses $M_{\rm 2}$, the orbital separations $a$ and the eccentricities. Since the eccentricity distribution has a minor effect on the BPS results (Hurley et al., 2002), we assume that all primordial binaries have circular orbits for simplicity. In our calculations, we take $M_{\rm 1}$ in the range of $5-100M_{\odot}$, $M_{\rm 2}$ in the range of $1-100M_{\odot}$, and $a$ in the range of $3-10^{4}R_{\odot}$. For the secondary stars, only the binaries with $M_{2}<M_{1}$ are included. We set each of the initial parameters $M_{1}$, $M_{2}$ and $a$ with the $n_{\chi}$ grid points of the parameter $\chi$ logarithmically spaced, so

Here $n_{\chi}$ is taken to be 100. If a specific binary $i$ evolves across a phase that is identified as a BH–CS system, then the binary contributes the BH–CS binary population with a rate

where $f_{\rm bin}$ is the binary fraction, $\rm SFR$ is the star formation rate, $M_{*}\sim 0.5M_{\odot}$ is the average mass for all stars, and $W_{\rm b}=\Phi(\ln M_{1})\varphi(\ln M_{2})\Psi(\ln a)\delta\ln M_{1}\delta\ln M_{2}\delta\ln a$ weights the contribution of the specific binary from the primordial binary with initial parameters of $M_{1}$, $M_{2}$ and $a$ (see more details in Hurley et al., 2002). We assume that all stars are initially in binaries, i.e. $f_{\rm bin}=1.0$. Since $\sim(60-90)\%$ of OB stars are observed as members of binary systems (Moe & Di Stefano, 2017), this assumption may lead to overestimate the population size of merging BH–CS binaries by a factor of less than 2. The primary masses are assumed to obey the initial mass function (Kroupa et al., 1993),

where $a_{1}=0.29056$ and $a_{2}=0.15571$ are the normalized parameters. Thus

The secondary masses are assumed to follow a flat distribution between 0 and $M_{1}$ (Kobulnicky & Fryer, 2007), then

The orbital separations are assumed to be uniformly distributed in the logarithm (Abt, 1983), thus

The normalization of this distribution yields $k=0.12328$.

We use the one-dimensional stellar evolution code MESA (version 10398, Paxton et al., 2011, 2013, 2015, 2018, 2019) to simulate the detailed evolution of a large grid of BH binaries with a nondegenerate companion star. The BH is treated as a point mass and its binary companion is initially a zero-age main sequence star. The evolutionary models are computed at solar metallicity ($Z_{\odot}=0.02$) and two sub-solar metallicities ($Z=$ 0.001 and 0.0001). The initial hydrogen mass fraction is assumed to be $X=1-Y-Z$, where $Y=0.24+2Z$ is the helium mass fraction (Tout et al., 1996). For the initial BH binaries, we take the BH masses distributed over a range of $3-20M_{\odot}$ in logarithmic steps of 0.3, the companion masses over a range of $1-100M_{\odot}$ in logarithmic steps of 0.05, and the orbital periods over a range of $1-10000$ days in logarithmic steps of 0.1.

In the MESA code, convective mixing is accounted for by using the mixing-length theory with a default mixing length parameter of $\alpha=2$. Following Brott et al. (2011), we include convective core-overshooting with an overshooting parameter of 0.335 pressure scale heights. Stellar winds are employed using mass-loss rate prescriptions similar to those in the BPS calculation. The wind mass-loss rates of Vink et al. (2001) are used for hot stars and the Nieuwenhuijzen & de Jager (1990) prescription for relatively cool stars with effective temperatures lower than $10^{4}$ K. For hydrogen-envelope stripped stars ($X<0.4$), we use the reduced rates of Hamann et al. (1995) for helium stars. We linearly interpolate the rates between different prescriptions to ensure a smooth transition as described by Brott et al. (2011). We adopt the scheme of Kolb & Ritter (1990) to calculate the mass-transfer rates $\dot{M}_{\rm tr}$ via RLO. Mass accretion onto the BH is limited by the Eddington rate $\dot{M}_{\rm edd}$, and the excess matter escapes from the binary system carrying away the specific orbital angular momentum of the BH (see also Shao & Li, 2020). We simulate the evolution by employing the default timestep options with $\rm mesh\_delta\_coeff=1.0$ and $\rm varcontrol\_target=10^{-4}$. Each binary evolution track is terminated if carbon is depleted in the companion’s core or the time steps exceed 30000.

Pavlovskii et al. (2017) showed that the mass-transferring BH binaries with a massive donor are likely to experience the expansion or the convective instability, and enter the CE evolution. The expansion instability occurs when the donor stars are experiencing a period of fast thermal-timescale expansion, while the convection instability takes place if the donor stars have developed a sufficiently deep convective envelope. Usually the former can be triggered in relatively close BH binaries, while the latter in relatively wide systems. It was suggested by Pavlovskii et al. (2017) that there exist the smallest radius $R_{\rm U}$ and the maximum radius $R_{\rm S}$ for which the convection and expansion instabilities can take place respectively. Mass transfer is stable if the donor radius is between $R_{\rm U}$ and $R_{\rm S}$. These modified CE criteria have been applied to the BPS study on the formation of double CS mergers (Olejak et al., 2021), but the limitation is that the calculated outcomes presented by Pavlovskii et al. (2017) are only from a small sample of binary evolution calculations.

For a BH binary with RLO mass transfer, the transferred material from the donor star forms an accretion disk around the BH. If the mass transfer proceeds at a very high rate, the matter will pile up around the BH and presumably form a bloated cloud engulfing a significant fraction of the accretion disk. Within a specific radius, photons can be trapped in the cloud. This is the so-called trapping radius that defined by Begelman (1979) as $R_{\rm trap}=(R_{\rm sch}/2)(\dot{M}_{\rm tr}/\dot{M}_{\rm edd})$, where $R_{\rm sch}$ is the Schwarzschild radius of the BH. We assume that the BH binary is engulfed in a CE if $R_{\rm trap}$ is larger than the RL radius $R_{L_{\rm BH}}$ of the BH, or the mass transfer rate exceeds a critical value, $\dot{M}_{\rm trap}=2\dot{M}_{\rm edd}R_{L_{\rm BH}}/R_{\rm sch}$ (King & Begelman, 1999; Belczynski et al., 2008). In other conditions, RLO via the $L_{1}$ point might not be rapid enough to remove the donor’s expanding envelope, so the donor star extends far beyond to reach the $L_{2}$ point. Mass loss through the $L_{2}$ point can take away a large amount of angular momentum from the binary system and lead to rapid binary orbit shrink, probably followed by the CE evolution (e.g. Ge et al., 2020; Misra et al., 2020). When dealing with the mass transfer stability in binaries with convective giant donors, Pavlovskii & Ivanova (2015) showed that the binaries can survive the mass transfer even after $L_{2}$ overflow without starting a CE phase. In our simulations, we assume that CE evolution takes place if either of the following conditions is met: (I) $\dot{M}_{\rm tr}>\dot{M}_{\rm trap}$; (II) $\dot{M}_{\rm tr}>0.02M_{\rm d}/P_{\rm orb}$ (Pavlovskii & Ivanova, 2015) and $R_{\rm d}>R_{L_{2}}$ (Ge et al., 2020). Here $M_{\rm d}$ is the donor mass, $P_{\rm orb}$ the orbital period of the BH binary, $R_{\rm d}$ the donor radius, and $R_{L_{2}}$ the volume-equivalent radius of the $L_{2}$ lobe. When CE starts, we terminate the calculation and record the relevant information for further analysis.

In Figure 2, we outline our calculated outcomes by showing the parameter space boundaries of the BH binaries between stable and unstable mass transfer. The panels from top to bottom correspond to the donor stars with different metallicities. The left panels show the boundaries for initial binaries in the orbital period $P_{\rm orb,i}$ vs. donor mass $M_{\rm d,i}$ diagrams, while the middle and right panels correspond to the boundaries for the binaries at the moment of starting RLO in the mass ratio $q$ (of the donor to the BH) vs. donor radius $R_{\rm d}$ and donor mass $M_{\rm d}$ vs. donor radius $R_{\rm d}$ diagrams, respectively. In each diagram, the coloured triangles represent the upper and lower boundaries for dynamically stable mass transfer with different initial BH masses. Our results are roughly coincident with those of Pavlovskii et al. (2017), showing that the binaries with the mass ratios up to $\sim 10$ can still proceed on a thermal timescale without entering CE evolution (see also Shao & Li, 2018; Marchant et al., 2021). In more detail, our simulations indicate that mass transfer in BH binaries is always stable if the mass ratio $q$ is smaller than the minimal value $q_{\rm min}$, i.e.

and always unstable if the mass ratio $q$ is larger than the maximal value $q_{\rm max}$, i.e.

where $M_{\rm BH}$ is the BH mass²²2We note that Equation (8) can also be applied to the binaries with an NS accretor, where the maximal mass ratio is $\sim 3-3.5$ (e.g., Kolb et al., 2000; Tauris et al., 2000; Podsiadlowski et al., 2002; Shao & Li, 2012; Misra et al., 2020).. For the binaries with mass ratio between $q_{\rm min}$ and $q_{\rm max}$, dynamically unstable mass transfer ensues if the donor radius is either less than $R_{\rm S}$, i.e.

or larger than $R_{\rm U}$, i.e.

Here all radii and masses are expressed in solar units. We roughly fit $R_{\rm S}$ and $R_{\rm U}$ as a function of $q$ (the orange dashed curve) and $M_{\rm d}$ (the orange dotted curve), respectively. When analysing our recorded data, we find that Equations (9) and (10) correspond to conditions (I) and (II), respectively. We also find that our obtained parameter spaces for stable mass transfer are not strongly dependent on stellar metallicities, except for the binaries with donors initially more massive than $\sim 40M_{\odot}$. At solar metallicity, mass transfer in long-period systems with very massive donors are always stable since the donor stars have experienced extensive wind mass loss prior to mass transfer (Klencki et al., 2021). As a consequence, we see in the $Z=0.02$ case that the green and blue filled triangles do not appear to show the upper boundaries for the binaries with $M_{\rm d,i}\gtrsim 40M_{\odot}$. Thus, we adopt Equations (7-10) to determinate whether the BH binaries enter the CE evolution in the BPS calculation. There are two exceptions: for the BH binaries with helium-star donors, we assume they can always avoid CE evolution (Tauris et al., 2015); for the BH binaries with WD donors, we assume mass transfer proceeds stably if the mass ratio of the WD to the BH is less than 0.628 (Hurley et al., 2002). When CE evolution is triggered, we allow donors that are crossing the Hertzsprung gap to survive the CE phase (the “optimistic” scenario of Belczynski et al., 2008).

Galactic BH$-$CS binaries are the high-mass analogues of the systems containing only NSs and/or WDs. The BH$-$WD and BH$-$NS systems may be discovered through observations of electromagnetic wave signals. On the other hand, close BH$-$CS binaries are likely to be observed by future space-based observatories due to the detection of GW signals. We will evaluate the population properties of close BH$-$CS binaries that can be identified by LISA.

The angle-averaged signal-to-noise ratio for a BH$-$CS system that can be detected over an observation time $T$ is given by (O’Leary et al., 2009)

where $n$ labels the harmonics at frequency $f_{n}\simeq nf_{\rm orb}$,

is the characteristic strain at the $n$th harmonic (Barack & Cutler, 2004), and $h_{N}(f_{n})$ is the characteristic LISA noise, including a contribution from unresolved Galactic binaries, for which we take from Robson et al. (2019). Here $T$ is taken to be $4\rm\,yr$ for the LISA mission duration, $f_{\rm orb}$ is the orbital frequency of the binary, $d_{\rm L}$ is the luminosity distance to the source, $G$ is the constant of gravity, and $c$ is the speed of light in vacuum. In Equation (12), $\dot{E}_{n}$ is the derivative of the energy radiated in GWs at frequency $f_{n}$, which to lowest order is given as (Peters & Mathews, 1963)

where $M_{\rm chirp}=(M_{\rm BH}M_{\rm CS})^{3/5}(M_{\rm BH}+M_{\rm CS})^{-1/5}$ is the chirp mass, $M_{\rm BH}$ and $M_{\rm CS}$ are the component masses of the BH$-$CS binary, and $g(n,e)$ is a function of the orbital eccentricity $e$ (from Peters & Mathews, 1963). To the leading quadrupole order, the term $\dot{f}_{n}$ is

where $F(e)=[1+(73/24)e^{2}+(37/96)e^{4}]/(1-e^{2})^{7/2}$. Note that a source is effectively less detectable the slower the GW frequency changes over the instrumental lifetime. This effect is taken into account by reducing the characteristic strain $h_{c,n}$ by a factor of the square root of $\min[1,\dot{f}_{n}(T/f_{n})]$ (see e.g., Kremer et al., 2019).

The peak frequency of GW emission for eccentric binaries is given by

where $n_{p}\simeq 2(1+e)^{1.1954}/(1-e^{2})^{1.5}$ (Wen, 2003). Obviously, $n_{p}=2$ for circular binaries. Considering that the majority of LISA-visible binaries have mild eccentricities, the GW power is sharply peaked at the peak frequency (Peters & Mathews, 1963). We follow Banerjee (2020) to simplify the signal-to-noise ratio as

with $n=n_{p}$ and $f_{n}=f_{\rm GW}$. We assume that a BH$-$CS binary can be identified by LISA if the signal-to-noise ratio is larger than 5 (e.g., Lamberts et al., 2018; Banerjee, 2020). To obtain the $d_{\rm L}$ distribution of the BH$-$CS binaries in the Milky Way, we assume that the binaries are uniformly distributed on a flat disc with a radius 15 kpc³³3Although this assumption is simple, it can roughly reflect the spatial distribution of detectable LISA binaries (see Figure 1 of Lau et al., 2020). and the Sun with respect to the Galactic Center has the distance of 8 kpc (Feast & Whitelock, 1997). When synthesizing the BH$-$CS binary population, we adopt continuous star formation at a rate ${\rm SFR}=3M_{\odot}\,\rm yr^{-1}$ (Smith et al., 1978; Diehl et al., 2006; Robitaille & Whitney, 2010) with solar metallicity over a period of 10 Gyr. In Table 1, we present the formation and merger rates of Galactic BH$-$CS binaries in our adopted three supernova models. Under all the models, we obtain that Galactic BH$-$WD, BH$-$NS and BH$-$BH binaries have the formation rates $\sim 11-95\rm\,Myr^{-1}$, $4-33\rm\,Myr^{-1}$ and $19-150\rm\,Myr^{-1}$, respectively, and the merger rates $\sim 0-3\rm\,Myr^{-1}$, $1-6\rm\,Myr^{-1}$ and $6-17\rm\,Myr^{-1}$, respectively. It can be estimated that the total number of Galactic BH$-$CS binaries formed over the past 10 Gyr is of the order $10^{5}-10^{6}$.

To estimate the merger rate density $\mathcal{R}_{\rm BHCS}$ of the BH$-$CS binaries in the local Universe, we have evolved a large number ($2\times 10^{7}$) of the primordial binaries for each model. The metallicity Z for every primordial binary is randomly taken in the logarithmic space between 0.0001 and 0.02. After simulations, we record relevant information for the BH$-$CS binaries that can merge within the Hubble time, including the CS types, the component masses and the delay time $t_{\rm delay}$ (from the formation of the primordial binaries to the merger of the BH$-$CS binaries). For each BH$-$CS merger, the look back time of the merger can be estimated as

where $t_{\rm lb}$ is the look back time for (binary) stars formed at redshift $z$,

where $H_{0}$, $\Omega_{\rm M}$ and $\Omega_{\Lambda}$ are the cosmological parameters for which we adopt the values from Planck Collaboration et al. (2016). Following Mapelli et al. (2017), we consider the BH$-$CS mergers with $t_{\rm merg}\geq 0$, excluding the systems that will merge in the future. And, we only include the mergers in the local Universe (defined as $z\leq 0.1$) using the condition of $t_{\rm merg}\leq t_{\rm lb}(z=0.1)$. We further divide the recorded binaries into 20 logarithmically spaced metallicity bins between $Z=0.0001-0.02$ for each model. The metallicity for stars at a given redshift is computed as $\log Z(z)/Z_{\odot}=-0.19z$ if $z\leq 1.5$ and $\log Z(z)/Z_{\odot}=-0.22z$ if $z>1.5$ (Rafelski et al., 2012). For $Z<0.0001$ or $Z>0.02$, we instead use the recorded information of the systems with $Z=0.0001$ or $Z=0.02$. In each metallicity bin, the parameter distribution of the primordial binaries has been normalized to unity. Similar to the treatment of Giacobbo & Mapelli (2018), we use the following analytic equation to calculate $\mathcal{R}_{\rm BHCS}$ as (see also Spera et al., 2019)

where $\mathcal{SFR}(z)$ is the cosmic star-formation rate density as a function of redshift for which we use the fitted formula given by Madau & Dickinson (2014),

In this paper, we have investigated the properties of merging BH$-$CS binary population in the Milky Way and the local Universe, based on binary evolution calculations with the BSE and MESA codes. Only the systems formed through isolated binary evolution are taken into account. Compared with previous works, the innovations in this work mainly lie in two aspects. (1) We have revised the criteria for the occurrence of CE evolution in the BH binaries with nondegenerate donors, by a large grid of detailed binary evolution simulations with the MESA code. As a consequence, we obtain the potential parameter space for dynamically (un)stable mass transfer, and incorporate into our BPS calculations. (2) We consider various mechanisms for the formation of NSs and BHs. The (non)existence of the mass gap between NSs and BHs is crucial to constrain the mechanism of supernova explosions. We adopt three models to deal with the compact remnant masses and natal kicks during supernova explosions, that is the rapid mechanism (Fryer et al., 2012), the delayed mechanism (Fryer et al., 2012) and the stochastic recipe (Mandel & Müller, 2020). The delayed and the stochastic models allow the formation of CSs within the mass gap, while the rapid model naturally leads to the mass gap between NSs and BHs. The fractions $f_{\rm MG}$ of merging systems with mass-gap BHs among the whole population for different types of BH$-$CS systems can be used as an indicator to examine relevant supernova mechanisms. In the rapid model, $f_{\rm MG}=0$ always holds.

We identify two formation channels for close BH$-$CS systems that appear as GW sources, depending on whether the mass transfer in the progenitor BH binaries is stable. We find that almost all close BH$-$WD binaries have experienced a CE phase, during which a mass-gap BH was engulfed by the envelope of a $\sim 6-10M_{\odot}$ supergiant (the WD’s progenitor). So for merging BH$-$WD systems, we expect that $f_{\rm MG}\sim 1$ in both the delayed and the stochastic models. For merging BH$-$NS and BH$-$BH systems, both the CE and stable mass transfer channels take effect. In both the delayed and stochastic models, we obtain that $f_{\rm MG}\sim 0.7$ for merging BH$-$NS binaries and $f_{\rm MG}\sim 0.3$ for merging BH$-$BH binaries. It seems that our $f_{\rm MG}$ predictions for merging BH$-$CS binary population cannot make a clear distinction between the delayed and the stochastic models. The LIGO/Virgo operation will provide a rapid growing sample of BH$-$BH merger events. For BH$-$BH mergers with total mass $\lesssim 30M_{\odot}$, the delayed model anticipates a more flat mass-ratio distribution between $\sim 0.4-0.9$ while the stochastic model favors large mass ratios distributed with a broad peak at $\sim 0.6-0.9$ (see Figure 12).

We estimate that there are totally dozens of BH$-$CS systems detectable by LISA in the Milky Way, and the merger rate density of BH$-$CS systems varies in the range of $\sim 60-200\rm\,Gpc^{-3}yr^{-1}$ in the local Universe. For merging BH$-$WD binaries, the delayed and stochastic models predict $\sim 2$ and $\sim 38$ systems may be observed by LISA in the Milky Way, respectively. In addition, we expect that local BH$-$WD merger rate is in the range of $\sim 7-59\rm\,Gpc^{-3}yr^{-1}$ for the delayed and stochastic models. Our calculations show that $\sim 2-14$ BH$-$NS systems are potential GW sources in the Milky Way and the local merger rate of BH$-$NS binaries is in the range of $\sim 10-80\rm\,Gpc^{-3}yr^{-1}$. Among all Galactic BH$-$BH systems, $\sim 12-26$ of them are expected to be the GW sources in the LISA frequency. Our calculated merger rate of BH$-$BH binaries in the local Universe varies in the range of $\sim 40-80\rm\,Gpc^{-3}yr^{-1}$. At last we remind that our results are still subject to many uncertainties, including the assumption of binary fraction ($f_{\rm b}=1$), the treatments of stellar and binary evolutionary processes (see e.g., Langer, 2012), the options of Galactic and cosmological parameters, and etc. It is obvious that $f_{\rm b}\sim 0.6-0.9$ from observations (Moe & Di Stefano, 2017) can lead to the decrease of the LISA-detectable numbers and the local rates for the merging BH–CS binary populations.