EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 651 (July 2021) A&A, 651 (2021) A100 Full HTML
Free Access
Issue
A&A
Volume 651, July 2021
Article Number A100
Number of page(s) 17
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202140520
Published online 27 July 2021

© ESO 2021

1. Introduction

The concept behind the formation of close binary systems, such as merging double compact objects (DCO) or X-ray binaries, via the common envelope (CE) phase, began to appear in the literature nearly 50 years ago (Paczynski 1976; van den Heuvel 1976). Based on this concept, one binary component enters the other binary companion’s envelope. The orbital energy is transferred to the envelope due to various drag forces, resulting in the binary orbit shrinking. Finally, the envelope either can be ejected from the system, leaving behind a newly formed close binary system, or the two stars can end up merging. One of the approaches to find the outcome of a CE event without doing detailed 3D simulations considers the energy budget of a CE event and is, hence, known as energy formalism. In its parameterized form, it has been introduced by Webbink (1984) and Livio & Soker (1988), and is known as “αCE-formalism”. Due to its simplicity, αCE-formalism is still widely used in population synthesis studies (e.g., Vigna-Gómez et al. 2020). This simplified form of energy formalism equates only two energies: the envelope’s binding energy (the energy required to eject the envelope to infinity) and the change in the orbital energy (the available energy source). This allows us to estimate the orbital separation after a CE event, assuming that envelope has been ejected:

G M don , i M don , env λ R don , lob = α CE ( G M don , f M comp , f 2 a f − G M don , i M comp , i 2 a i ) . $$ \begin{aligned} \frac{G M_{\rm don,i} M_{\rm don,env}}{\lambda R_{\rm don,lob}} = \alpha _{\rm CE} \left(\frac{G M_{\rm don,f} M_{\rm comp,f}}{2 a_{\rm f}} -\frac{G M_{\rm don,i} M_{\rm comp,i}}{2 a_{\rm i}} \right). \end{aligned} $$(1)

Here, ai/f is the initial or final orbital separation, Mdon, i/f – initial or final donor mass, Mcomp, i/f – initial or final companion mass, respectively, G – gravitational constant, Mdon, env – mass of the donor envelope, Rdon, lob – Roche lobe radius of the donor at the onset of Roche-lobe Overflow (RLOF), i and f – initial and final values of mass and separation, and λ – measure of the donor central concentration (de Kool 1990; Dewi & Tauris 2000; Xu & Li 2010).

Equation (1) postulates that the transfer of binary system orbital energy into the energy of the envelope takes place with some efficiency, αCE, which, unless other energy sources are present, can not be more than one. Simulations and observations of CE phase indicates that the value of parameter is typically αCE <  0.6−1.0 (Zuo & Li 2014; Nandez & Ivanova 2016; Clayton et al. 2017; Iaconi & De Marco 2019). On the other hand, if the simulations of a CE event are performed while including more physical processes, for example, accretion, the effective value of αCE can be as high as 5 (e.g., see Fragos et al. 2019). A similar effect of increasing the apparent CE efficiency to more than one can be produced via exotic nucleosynthesis and jets (Podsiadlowski et al. 2010; Shiber et al. 2019; Zevin et al. 2021; Grichener & Soker 2021). Some population synthesis studies have already adopted such high values (e.g., see Santoliquido et al. 2020).

The CE efficiency parameter αCE and the binding energy parameter λ are often coupled in population synthesis studies or in observation analyses when energy formalism is applied. At the same time, each of them is subject to many uncertainties. For example, λ was introduced to relate the “true” binding energy of the envelope (as obtainable from detailed stellar models) to its simple parameterized form. However, the true binding energy of the donor continues to be a subject of discussion (e.g., see Sect. 3 in Ivanova et al. 2020). The CE efficiency parameter is expected to depend on the system specifics and which physical processes of energy creation or energy loss took place (e.g., on the system’s mass ratio, the evolutionary stage of the donor, or the nature of the companion). Values for αCE, derived from observations, seem to be systematically lower for asymptotic giant branch (AGB) donors than for red giant branch (RGB) donors (Iaconi & De Marco 2019), which have different internal structures and envelope binding energy. The AGB donor envelopes are considered to be less tightly bound to the core than those of a RGB donors (Paczyński & Ziółkowski 1968; Han et al. 1994). However, even taking into account this AGB stars feature does not sufficiently help in successive envelope ejection (Sand et al. 2020). Numerical simulations encounter difficulties in successful envelope ejection unless some other energetic process (except orbital energy release) is included (e.g., Passy et al. 2012; Iaconi et al. 2017; Sand et al. 2020). In addition, recent studies by Klencki et al. (2021) have confirmed that even with the most favorable assumptions, a successful CE ejection in BH binaries is only possible if the donor is a massive convective-envelope giant. On the other hand, massive stars (BH progenitors) may be a subject of extensive mass loss through enhanced winds before they reach the RSG stage, which may even cause the spontaneous envelope loss (e.g., Vanbeveren 1991; Vanbeveren et al. 1998; Eggleton 2002).

The energy formalism in its parameterized form is a convenient way to predict CE outcomes, but, as it is argued now, it is not necessarily a well-founded method (Ivanova et al. 2020). While, at the moment, there is no comprehensive understanding of the CE evolution in all the cases, it is understood that the CE event could be preceded by stable mass transfer (MT) on different timescales. The criteria for the occurrence of CE phase are still under development as well.

New stellar mass loss models (Ge et al. 2010, 2015, 2020a; Ge et al. 2020b) and detailed simulations for close, mass-exchanging binaries (Pavlovskii et al. 2017; Misra et al. 2020) have shown that RLOF may be stable over a much wider parameter space than previously thought. The same studies indicate that RLOF stability depends not only on the system mass ratio and the envelope type (convective or radiative), but also, for instance, on the metallicity, stellar type, or radius of the donor star. Those results can be confirmed via comparisons of the theoretical models with the observed systems (Cherepashchuk et al. 2019; Leiner & Geller 2021). A summary of the recent progress on RLOF stability has been summarized in Sect. 2.2 of Klencki et al. (2021). Unfortunately, due to a scarcity of observations for massive stellar systems during the ongoing RLOF phase and the high calculation costs, the current CE study (observations and simulations) usually refers to low-mass binary systems (Nandez et al. 2015; Nandez & Ivanova 2016; Jones 2020), which are not progenitors of BH-BH, BH-NS, or NS-NS binaries.

The CE phase is a key element in setting the formation of DCOs in the classical isolated binary evolution channel and therefore understanding of CE is crucial in studies of origin of merging DCOs. Recently more and more signals form BH-BH, BH-NS and NS-NS mergers have been detected by LIGO/Virgo instruments (Abbott et al. 2019a; The LIGO Scientific Collaboration 2021) during the O1, O2 and O3 runs with reported parameters of systems such as the masses, spins and redshifts. It is still unknown what fraction of gravitational wave (GW) signal mergers formed through isolated binary evolution in the field (Bond & Carr 1984; Tutukov & Yungelson 1993; Lipunov et al. 1997; Voss & Tauris 2003; Belczynski et al. 2010; Dominik et al. 2012; Kinugawa et al. 2014; Mennekens & Vanbeveren 2014; Hartwig et al. 2016; Spera et al. 2016; Belczynski et al. 2016b; Eldridge & Stanway 2016; Woosley 2016; Stevenson et al. 2017; Kruckow et al. 2018; Hainich et al. 2018; Marchant et al. 2019; Spera et al. 2019; Bavera et al. 2020), the dense stellar system dynamical channel (Miller & Hamilton 2002b,a; Portegies Zwart et al. 2004; Gültekin et al. 2004, 2006; O’Leary et al. 2007; Sadowski et al. 2008; Downing et al. 2010; Antonini & Perets 2012; Benacquista & Downing 2013; Bae et al. 2014; Chatterjee et al. 2017; Mapelli 2016; Hurley et al. 2016; Rodriguez et al. 2016; VanLandingham et al. 2016; Askar et al. 2017; Morawski et al. 2018; Banerjee 2018; Di Carlo et al. 2019; Zevin et al. 2019; Arca-Sedda & Capuzzo-Dolcetta 2019; Rodriguez et al. 2018; Perna et al. 2019; Kremer et al. 2020); isolated multiple (triple, quadruple) systems (Antonini et al. 2017; Silsbee & Tremaine 2017; Arca-Sedda et al. 2021; Liu & Lai 2018; Fragione & Kocsis 2019), mergers of binaries in galactic nuclei (Antonini & Perets 2012; Hamers et al. 2018; Hoang et al. 2018; Fragione et al. 2019); the chemically homogeneous evolution channel consisting of rapidly spinning stars in isolated binaries (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016; du Buisson et al. 2020), or Population III origin DCO binary mergers (Bond & Carr 1984; Kinugawa et al. 2014; Tanikawa et al. 2021).

Over the years, various groups have developed their population synthesis codes to try and put better constrains on astrophysical processes by comparing theoretical model results with the known Galactic and extragalactic compact object population (e.g., Portegies Zwart & Verbunt 1996; Hurley et al. 2002; Kruckow et al. 2018; Belczynski et al. 2020; Han et al. 2020). Several uncertainties commonly encountered in population synthesis studies strongly influence formation of DCOs and therefore, are subjects of active research. Examples of uncertain processes and parameters that are crucial for DCO mergers evolution are: the metallicity-specific star formation rate density (Chruślińska et al. 2020; Santoliquido et al. 2021; Broekgaarden et al. 2021), NS and BH natal kicks (Mandel et al. 2021), or MT during stable or unstable RLOF (Vinciguerra et al. 2020; Howitt et al. 2020; Bavera et al. 2020).

In this paper, we study how the application of the most recent developments in the field of RLOF stability to the population synthesis affects DCOs formation. In Sect. 2 we describe the general method and the input physics implemented in the current version of StarTrack population synthesis code. In Sect. 3 we introduce the revised CE development criteria and modified stable RLOF treatment examined in this paper. In Sect. 4 we present the results of our simulations: DCO local merger rate density, BH-BH and BH-NS mass ratio distributions, and BH-BH mass distributions for three tested models. Section 5 is a description of evolutionary scenarios leading to the BH-BH and BH-NS mergers formation in three tested models. It also includes three diagrams with an example of systems evolution. Section 6 contains a brief discussion of our results together with the conclusions.

2. Method

In our simulation of the formation and mergers of DCO systems in the local Universe, we used the updated StarTrack population synthesis code, developed over a number of years (Belczynski et al. 2002, 2008). The code allows us to simulate the evolution of a single star as well as of a binary star system for a wide range of initial conditions and physical parameters. The currently used version of implemented physics, adopted star formation history (SFH), and metallicity of the Universe is described in Belczynski et al. (2020), with two recent modifications, both explained in Sect. 2 of Olejak et al. (2020). We adopted three broken power-law initial mass function (IMF) Kroupa et al. (1993, 2002), weak pulsation pair-instability supernovae (PPSN), and pair-instability supernovae (PSN) (Woosley 2017; Belczynski et al. 2016a). We applied procedures for accretion onto a compact object during stable RLOF and from stellar winds, based on the analytic approximations described in King et al. (2001) and Mondal et al. (2020). For non-degenerate accretors, we adopted a 50% non-conservative RLOF (Meurs & van den Heuvel 1989; Vinciguerra et al. 2020) with a fraction of the lost donor mass accreted onto the companion (fa = 0.5), and the rest of the mass (1 − fa), leaving the system together with part of the donor and orbital angular momentum (see Sect. 3.4 of Belczynski et al. 2008). We use 5% Bondi-Hoyle rate accretion onto the compact object during the CE phase (Ricker & Taam 2008; MacLeod & Ramirez-Ruiz 2015; MacLeod et al. 2017). The procedure is based on Eqs. (5.3)–(5.7) of (Bethe & Brown 1998) and Eqs. (A.1)–(A.10) of Belczynski et al. (2002), and has been recently summarized in Appendix B of Olejak et al. (2020). For stellar winds we use formulas based on theoretical predictions of radiation driven mass loss (Vink et al. 2001) with the inclusion of Luminous Blue Variable mass loss (Belczynski et al. 2010). All tested models have our standard physical values for the envelope ejection efficiency of αCE = 1.0 and Maxwellian distribution natal kicks with σ = 265 km s−1 (Hobbs et al. 2005) lowered by fallback (Fryer et al. 2012). We adopted solar metallicity Z = 0.02, in agreement with Pavlovskii et al. (2017).

We assume that systems with Hertzsprung gap (HG) donor star merge during CE phase (Belczynski et al. 2007). In the StarTrack code, the HG phase begins after leaving the main sequence – for both less and more massive stars, it is the period of intense star expansion (since, during the main sequence, stellar radius usually does not increase more than by a factor of few). Therefore, during the HG phase, stars often initiate TTMT and CE. At the onset of RLOF, such donors are often only partially expanded post-main sequence stars. Therefore, it is not well known whether such objects have an already well-separated core and envelope structure. In Fig. 1, we present the Hertzsprung-Russell evolution diagram for massive single stars. In the top panel of the figure, we mark that part of the evolution when stars are expected to have radiative or convective envelopes. In the bottom panel, we mark the part of the evolution where stars are defined in StarTrack as main sequence+HG or core helium burning.

thumbnail Fig. 1.

Hertzsprung-Russell diagram for massive stars with ZAMS masses: 20 M⊙, 32 M⊙, 50 M⊙, 80 M⊙, and 100 M⊙ at metallicity Z = 0.007; evolution of single star since main sequence to core helium burning. In the top panel, we marked stars cooler than log(Teff/K) < 3.73 with the red line and with the blue line, those stars with log(Teff/K) > 3.73. It is expected that the red line represents stars having convective envelopes, while the blue line is for radiative envelopes (Ivanova & Taam 2004; Belczynski et al. 2008). In the bottom panel, we marked with the blue line, the star main sequence and Hertzsprung gap phases and with the red line, the core helium burning phase. Evolutionary stages are defined as in formulas by Hurley et al. (2002) and Belczynski et al. (2008). With the grey lines, we marked the borders corresponding to stellar radii: 100 R⊙, 500 R⊙, and 1000 R⊙.

We tested a delayed supernovae (SN) engine (Fryer et al. 2012; Belczynski et al. 2012), which affect the birth mass of NSs and BHs, allowing for the formation of the compact objects within the first mass gap (∼2−5 M⊙). We assume that maximum mass of the NS is 2.5 M⊙ (Horvath et al. 2020) so more massive compact objects are BHs.

In this work, we present our results for three different evolutionary models: a model with standard Startrack treatment of RLOF: M380.B, along with the other two, which include the revised RLOF treatment: M480.B and M481.B. In these names, “M” stands for model, then the number is increasing with time to mark subsequent physical models. The ending “B” refers to submodel B, which is often used in previous Startrack works, such as in Belczynski et al. (2020). In submodel B, all HG donor systems merge during the CE phase. The models are listed in Table 1.

Table 1.

List of tested evolutionary models.

3. RLOF

Two new types of instabilities were identified in RLOF binaries with massive donors, if such donors were allowed to evolve while they are exceeding their Roche lobe – the expansion and the convective instability (Pavlovskii et al. 2017). Expansion instability happens if an RLOF donor experiences a period of fast thermal-timescale expansion after its main sequence. This expansion may lead to the development of dynamical instability in a very short time (a few thousand years) after the start of thermal time-scale mass transfer (TTMT). The second type of instability, convection instability, is associated with developing a sufficiently deep convective envelope.

In this work, we implement, in our population synthesis code, the revised development criteria for the occurrence of the two instabilities and we examine how they influence the formation of DCO mergers. We base our criteria for checking whether any of the instabilities occur in a given donor on the numerical results of Pavlovskii et al. (2017), who report the boundaries between the instabilities using the radii of the donors of different masses and metallicities at the onset of RLOF. Following Pavlovskii et al. (2017), we defined two radii: RU, the smallest radius for which the convection instability occurs, and RS, the maximum radius when the expansion instability can take place.

In our approximation, the values of RU and RS are the averages of the ranges obtained and given in Table 1 of Pavlovskii et al. (2017). The specific values of RU and RS used in this work are given in Table 2. We note that simulations by Pavlovskii et al. (2017) were performed for the systems with BH companion. In our work, we extend the results and treat the same way binary systems with any other companion types (e.g., main-sequence stars, neutron stars, etc.). If the radius, Rdon, of a given donor at the onset of RLOF is found within the range RS–RU for Z ≤ 0.5 Z⊙ or below Rdon <  RU for Z >  0.5 Z⊙, then we assume a stable RLOF; otherwise, CE evolution is applied. We note that stars of a higher metallicity were not found to undergo such an expansion instability (Pavlovskii et al. 2017).

Table 2.

Boundary radii values between the stable and unstable RLOF, values based on data from Pavlovskii et al. (2017).

3.1. CE development criteria

In the revised treatment, we require the following four conditions to be met simultaneously for CE development:

The condition based on the donor type. Revised RLOF instability criteria are applied if, during the RLOF, the donor is an H-rich envelope giant of one of the following types: HG star, first giant branch, core helium burning, early asymptotic giant branch, or thermally pulsing asymptotic giant branch. However, for HG donors, we assume that CE always leads to a merger, and we halt this system binary evolution (see Sect. 2).

The condition based on the donor mass. Revised RLOF instability criteria are applied only to the systems with initially massive donor stars with MZAMS, don >  18 M⊙.

The condition based on the system mass ratio. The ratio of companion masses, Mcomp, to donor masses Mdon (see Table 2) represents the border between regimes of always stable and possibly unstable RLOF. Therefore, from Table 2 we can obtain the limit value of mass ratio, qCE = Mcomp/Mdon, for the binary to enter a potentially unstable RLOF regime. For metallicities Z ≤ 0.5 Z⊙ and for different donor mass, Mdon ranges the qCE values are following:

q CE = 0.36 for M don ∈ ( 18 , 60 ) M ⊙ q CE = 0.21 for M don ∈ [ 60 , 80 ) M ⊙ q CE = 0.19 for M don ≥ 80 M ⊙ , $$ \begin{aligned}&{q_{\rm CE}} = 0.36\ \mathrm{for }\ M_{\rm {don}} \in (18,60) M_{\odot } \nonumber \\&{q_{\rm CE}} = 0.21\ \mathrm{for }\ M_{\rm {don}} \in [60,80) M_{\odot } \nonumber \\&{q_{\rm CE}} = 0.19\ \mathrm{for }\ M_{\rm {don}} \ge 80 M_{\odot }, \end{aligned} $$(2)

whereas for metallicities Z >  0.5 Z⊙

q CE = 0.36 for M don ∈ ( 18 , 60 ) M ⊙ q CE = 0.29 for M don ∈ [ 60 , 80 ) M ⊙ q CE = 0.19 for M don ≥ 80 M ⊙ . $$ \begin{aligned}&{q_{\rm CE}} = 0.36\ \mathrm{for}\ M_{\mathrm{don}} \in (18,60) M_{\odot } \nonumber \\&{q_{\rm CE}} = 0.29\ \mathrm{for}\ M_{\mathrm{don}} \in [60,80) M_{\odot } \nonumber \\&{q_{\rm CE}} = 0.19\ \mathrm{for}\ M_{\mathrm{don}} \ge 80 M_{\odot } . \end{aligned} $$(3)

If the mass ratio of the binary system is greater than qCE, we assume stable RLOF.

The condition based on the donor radius. Finally, if all four conditions are met, we use stability diagrams presented in Figs. 2 and 3 to decide between stable RLOF and CE development. The diagrams pertain to two metallicities corresponding to the simulations performed by Pavlovskii et al. (2017): Z = 0.1 Z⊙ and Z = 1.0 Z⊙, respectively. We marked the values of RU and RS with the corresponding donor masses taken from the Table 2 and fit functions to the simulation points. We extend our model for wider metallicity ranges and follow different procedures for metallicities smaller than 0.5 Z⊙ and greater than 0.5 Z⊙.

thumbnail Fig. 2.

Revised CE development criteria for stars with 10% of solar metallicity. On horizontal axis is the mass of the donor star. On the vertical axis is the radius of the donor. The circle and square points mark values of RS and RU, respectively (see Table 2). We fit lines to these data (see Eqs. (4)–(6)) to show regions of stable RLOF and unstable RLOF (CE). The color coded bar on the right denotes companion mass.

thumbnail Fig. 3.

Revised CE development criteria for stars with 100% of solar metallicity. Notation is the same as on Fig. 2.

For Z ≤ 0.5 Z⊙, we used simulation data for Z = 0.1 Z⊙ and we fit a second-degree polynomial to RS points with the following coefficients:

R S = − 0.29 M don 2 + 30.3 M don − 498 , $$ \begin{aligned} {R_{\rm {S}} = -0.29 M_{\rm {don}}^2+30.3 M_{\rm {don}}-498}, \end{aligned} $$(4)

and a straight line to RU:

R U = 26.3 M don + 262 . $$ \begin{aligned} {R_{\rm {U}} = 26.3 M_{\rm {don}}+262.} \end{aligned} $$(5)

If a donor during ongoing RLOF has a radius Rdon < RS or Rdon > RU, RLOF is unstable (CE).

For Z >  0.5 Z⊙ we used simulation data for Z = 1.0 Z⊙. In this case, we had only RU points and we fit a straight line with the equation:

R U = 62.3 M don − 515 . $$ \begin{aligned} {R_{\rm U} = 62.3 M_{\rm don}-515.} \end{aligned} $$(6)

If a donor radius during ongoing RLOF is Rdon > RU, RLOF is unstable (CE).

Mdon in Eqs. (4)–(6) stands for the mass of the donor at a given RLOF time step, such that only 1% of donor mass is transferred (see Sect. 5.1, Belczynski et al. 2008).

In detailed simulations, RLOF is always unstable if q <  0.123 (see also Table 2). We therefore adopt the premise that CE events always start in binaries where at the onset of RLOF the mass ratio is

M comp M don < q crit = 0.125 . $$ \begin{aligned} \frac{M_{\rm comp}}{M_{\rm don}} < q_{\rm crit}=0.125. \end{aligned} $$(7)

Unfortunately, the He giant donors do not have envelopes similar to hydrogen stars to make a straightforward connection and use revised RLOF stability diagrams. Therefore, in this work, for He stars – as well as cases where the donor is a main sequence star (H-rich; core H-burning and He-rich; core He-burning main sequence) or a giant of MZAMS, don <  18 M⊙ – we follow the standard StarTrack CE development criteria and stable RLOF treatment (Belczynski et al. 2008).

3.2. Stable RLOF

In this study, we also test a modification of the condition for switching from TTMT to nuclear-timescale stable MT in the StarTrack code. So far, once the TTMT began, the issued MT was continued until the donor star’s radius Rdon decreased below 1.1 of its Roche lobe radius Rlobe. Then the type of MT was selected based on the comparison of MT timescale τeq (Eq. (46) from Belczynski et al. 2008) and thermal timescale τthermal. For τthermal, we use the approximation from Kalogera & Webbink (1996):

τ thermal ≈ 30 ( M don M ⊙ ) 2 ( R don R ⊙ ) − 1 ( L don L ⊙ ) − 1 Myr . $$ \begin{aligned} \tau _{\rm thermal} \approx 30 \left(\frac{M_{\rm don}}{{\,M}_{\odot }}\right)^2 \left(\frac{R_{\rm don}}{{\,R}_{\odot }}\right)^{-1}\left(\frac{L_{\rm don}}{{\,L}_{\odot }}\right)^{-1}\mathrm{Myr} . \end{aligned} $$(8)

If τeq >  τthermal, it was assumed that MT is driven on nuclear timescale. If τeq ≤ τthermal, it was assumed that the MT proceeds on a thermal timescale.

In one of the tested models of this work, we modified this condition by allowing for the switch from TTMT to nuclear timescale stable MT that is always based on the timescale comparison in the given RLOF time step (not only if Rdon <  1.1Rlobe). We note that this modification is applied to all binaries during their RLOF, regardless of component masses or types (as opposed to the new CE development criteria described in Sect. 3.1).

The evolutionary formulas we employ to evaluate star properties (e.g., radii, mass, luminosity) based on Hurley et al. (2000) are not applicable for thermal-timescale changes of the star outside of thermal-timescale evolution that is intrinsic to normal single stars (e.g., post-main-sequence expansion). Therefore, we cannot properly evaluate the donor’s properties during ongoing TTMT, as then the donor is out of its thermal equilibrium. We obtain TTMT rate using the mass, radius, and luminosity that the donor had when the TTMT had just started. Hence, the donor is still in thermal equilibrium. Therefore, we assume that the donor’s properties should still be well approximated by our parameterized equations. This obtained TTMT rate is kept constant for as long as the TTMT proceeds. We have to resort to this simplified method of obtaining the MT rate since, as the TTMT proceeds, our evolutionary formulas generate unreliable values for the star’s radius and luminosity.

However, we go on to use the onset RLOF parameters to construct our switches. In the first (original StarTrack) switch from TTMT to nuclear timescale MT, we wait until the donor’s radius (as given by the evolutionary formula) becomes comparable to the size of the donor’s Roche lobe. Then we assume that the star is back (or close) to its thermal equilibrium. Then we use the new donor’s properties to decide whether to continue TTMT (with a new MT rate based on the current donor mass, radius, and luminosity) or to switch to nuclear timescale MT.

In the modified approach that we test in this work, at each time step of the ongoing TTMT, we use the evolutionary formulae and calculate the current thermal and nuclear timescales. Suppose we get that the updated nuclear timescale is longer than the updated thermal timescale and the donor radius is still larger than its Roche lobe radius. In that case, we assume that the donor has regained its thermal equilibrium and we proceed with the MT on the new nuclear timescale. This usually happens earlier in the ongoing TTMT than meeting the original condition, Rdon <  1.1Rlobe. With this new additional condition, a numerical artifact is plausible: the donor star’s radius may significantly exceed its Roche lobe when we change MT to nuclear timescale. In this case we apply a “safety” condition, such that if:

R don > 2.0 R lobe , $$ \begin{aligned} R_{\rm don} > 2.0 R_{\rm lobe}, \end{aligned} $$(9)

we proceed with CE evolution for giant donors, or assume merger for main sequence donors (both H-core burning H-rich and He-core burning He-rich main sequence stars).

Using either the original switch or the modified one is not perfect. Since we cannot be sure which transition condition from TTMT to nuclear timescale MT would correspond better to detailed MT calculations, we will present models with both treatments.

3.3. Other cases of unstable RLOF

There are other situations when RLOF is unstable, leading to a CE phase. We take into account the following two of the known scenarios: (i) in case of the Darwin instability, due to an extreme mass ratio (Lai et al. 1993); (ii) if the accretion flow is so strong that photons are trapped in it, building up an envelope around a compact accretor in excess of its Roche lobe radius (Begelman 1979). For these, we use standard StarTrack procedures (Belczynski et al. 2008).

4. Results

4.1. NS-NS, BH-NS, and BH-BH local merger rate density

We present local merger rate densities (z ∼ 0) for different types of DCO systems corresponding to the tested physical models with standard (M380.B) and revised RLOF treatment (M480.B, M481.B) listed in the Table 3. Our rates are placed together with values recently estimated by LIGO/Virgo (90% credible limits): 23 . 9 − 8.6 + 14.9 \,Gpc − 3 \,yr − 1 $ 23.9^{+14.9}_{-8.6} {\text{\,Gpc}^{-3} \text{\,yr}^{-1}} $ for BH-BH mergers, 320 − 240 + 490 $ {320^{+490}_{-240}} $ Gpc−3 yr−1 for NS-NS mergers (The LIGO Scientific Collaboration 2021) and previously given range for BH-NS mergers: 0−610 Gpc−3 yr−1 (Abbott et al. 2019a).

Table 3.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different types of DCO systems.

The effect of revised CE development criteria (separate from other tested changes) on the merger rates is visible by comparing results of model M380.B and M480.B. It seems to be rather counterintuitive as the rate for BH-BH and BH-NS mergers increased under more restricted conditions for CE initiation. In the standard model (M380.B), the local merger rate densities were for BH-BH: 62 Gpc−3 yr−1 and for BH-NS: 13 Gpc−3 yr−1 while in M480.B they slightly increased to 88 Gpc−3 yr−1 and 16 Gpc−3 yr−1, respectively. This is caused by the fact that in the revised treatment during the RLOF some of the systems with HG star donor instead of entering the unstable RLOF regime (which always lead to the merger, Sect. 2), go through the period of TTMT which allows for further binary evolution and potential formation of BH-BH systems. On the other hand, systems initiate CE phase much less often and as the consequence we obtain new, dominant evolutionary scenario leading to formation of BH-BH mergers without any CE phase. In other words, the early (HG donors) development of CE during TTMT eliminates many progenitors of BH-BH systems from dominant formation channel of model M480.B and the restricted (in terms of mass and mass ratio) development of CE reduces number of BH-BH mergers from the dominant formation channel in model M380.B. For details about the new evolutionary scenario see Sect. 5. The rate for NS-NS mergers (model M480.B), ∼150 Gpc−3 yr−1, is the same as in the standard model (M380.B), as the revised CE development criteria were adopted only to the systems with initially massive donors (MZAMS >  18 M⊙, see Sect. 3.1) which are mainly BH progenitors. Therefore NS-NS systems formation was not significantly influenced and the merger rate is within LIGO/Virgo range as it used to be for the standard approach. Rates for BH-NS systems for both M380.B and M480.B models are also consistent with the given LIGO/Virgo range, which is wide, due to the limited detection data for BH-NS mergers (The LIGO Scientific Collaboration 2021). BH-BH merger rates for the same two models are about 1.5–2.0 times larger than the upper limit given by LIGO/Virgo.

Model M481.B includes both implemented changes to our standard approach (M380.B): revised CE development criteria and modified condition for the switch between TTMT and nuclear timescale stable MT. As the second change was implemented in all systems (regardless components masses or types), it affected also NS-NS merger rate (for M481.B equal ∼320 Gpc−3 yr−1), increasing it about twofold in comparison with the M380.B and M480.B models. The additional modification implemented in model M481.B leads, however, to a decrease in the rate for BH-BH and BH-NS mergers (18 Gpc−3 yr−1 and 4 Gpc−3 yr−1 respectively) by a factor of ∼4−5 comparing to the model M480.B. The effect described in previous paragraph – leading to the increase in the number of the surviviving BH binaries progenitors due to the entry into the stable RLOF instead of early CE with HG donors – does not apply here. Due to the modification of condition for switch between MT types, TTMT stops faster than in the standard approach. At the same time the radius of the donor, derived from the formula for stars in equilibrium (see Sect. 3.2) is large (Rdon >  2.0Rlobe), leading to early CE initiation (and merger of all HG donor systems). The new dominant BH-BH merger formation channel emerges in this model and it is different from the dominant BH-BH formation channels in models M380.B and M480.B. This channel includes TTMT that develops into late CE with a core helium burning donor (for details see Sect. 5). Merger rates for all types of DCOs in this (M481.B) model are in good agreement with recent LIGO/Virgo estimates.

4.2. BH-NS and BH-BH mass ratio distribution

Until the O3 data became available, all ten detected O1/O2 BH-BH mergers published by LIGO/Virgo team were consistent with being equal-mass mergers (Abbott et al. 2019b,c; Fishbach & Holz 2020). The latest detections, however, have indicated that the distribution of BH-BH mergers mass ratio may be more complex and wide-ranging. Recently, two detections of DCO mergers with highly asymmetric masses were announced (The LIGO Scientific Collaboration & the Virgo Collaboration 2020; Abbott et al. 2020a). The events GW190412 (∼30 M⊙ BH and ∼8 M⊙ BH) and GW190814 (∼23 M⊙ BH and ∼2.6 M⊙ NS or BH) were caused by DCO mergers in which one of the objects was less massive than the second by a factor of 0 . 28 − 0.07 + 0.12 $ 0.28_{-0.07}^{+0.12} $ and 0 . 112 − 0.009 + 0.008 $ 0.112_{-0.009}^{+0.008} $, respectively. We check how revised CE development criteria and stable RLOF treatment influence the mass ratio distribution of BH-BH, and BH-NS mergers, and if we are still able to reconstruct asymmetric mass BH-BH and BH-NS mergers as in previous models (Olejak et al. 2020; Drozda et al. 2020). In our study, we define the mass ratio of merging components, q, as the ratio of the less massive compact object to the more massive compact object, in agreement with LIGO/Virgo definition. In Figs. 4 and 5 we present the distributions of mass ratio of BH-BH and BH-NS mergers (z ∼ 0) for two models with revised RLOF treatment (red line): M480.B and M481.B, respectively.

thumbnail Fig. 4.

Distribution of mass ratio in BH-BH and BH-NS mergers (z ∼ 0). Top panel: BH-BH mergers, on the bottom panel BH-NS mergers. Blue line – results for standard CE development criteria with delayed SN engine (M380.B). Red line – results for revised CE development criteria (M480.B).

thumbnail Fig. 5.

Distribution of mass ratio in BH-BH and BH-NS mergers (z ∼ 0). On the top panel BH-BH mergers, on the bottom panel BH-NS mergers. Blue line – results for standard CE development criteria with delayed SN engine (M380.B). Red line – results for revised CE development criteria and revised stable RLOF treatment (M481.B).

We calculated the percent of asymmetric mass BH-BH mergers with q <  0.4, which has been constrained after GW190412 detection to constitute ≳10% of the overall BH-BH merger population (The LIGO Scientific Collaboration & the Virgo Collaboration 2020). In the standard StarTrack models with rapid (M230.B, Olejak et al. 2020) and delayed (M380.B) SN engine (Fryer et al. 2012), those fractions are 9% and ∼30%, respectively. We find that our tested revised CE development criteria reduced the fraction of unequal BH-BH mergers in comparison to model M380.B. Percent of BH-BH mergers with q <  0.4 for model M480.B constitute ∼2% while for M481.B it is ∼6%. The main formation channel for the formation of unequal mass BH-BH binaries in standard approach (see Table 1 of Olejak et al. 2020) is missing for the revised CE treatment (see Table 4, Sect. 5). Systems that would produce unequal mass BH-BH mergers, in the revised treatment instead of CE, enter the stable RLOF regime and finally form wide BH-BH binaries, which do not merge in Hubble time. The lowest BH-BH merger mass ratio achieved in tested models is q ∼ 0.06 for the M480.B and q ∼ 0.1 for M481.B.

In distribution of BH-BH mergers for model M480.B there is a large peak in the range of q ≈ 0.4−0.6, which is not present in distributions for other models. Vast majority (94%) of BH-BH mergers in M480.B forms via evolutionary scenario without any CE phase (see Table 4 and Fig. 9). Binary separation, instead during CE phase as in the standard Startrack approach, is reduced as the donor mass is ejected from the system together with the orbital angular momentum during stable RLOF. The described mechanism of orbital shrinkage is effective only for unequal mass binary components, when the donor mass at the RLOF onset is significantly larger than the companion’s mass (van den Heuvel et al. 2017; Marchant et al. 2021). In our simulations such systems usually end their evolution as BH-BH binaries with mass ratio close to q = 0.5 (Fig. 9).

Table 4.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different BH-BH formation scenarios.

Mass ratio distribution for BH-NS mergers is strongly dominated with low mass ratio systems, what has been already reported and described for example by Drozda et al. (2020). For both models M480.B and M481.B, ∼90% of the BH-NS mergers have q <  0.35. However, the lowest achieved mass ratios for BH-NS merger (∼0.1) is not that extreme as it used to be in standard model M380.B (0.03). Similarly as in case of the BH-BH system, the main formation channel leading to the formation of extreme mass ratios (see Fig. 6 in Drozda et al. 2020) is reduced as during the first RLOF instead of CE systems goes through period of stable RLOF.

4.3. BH-BH mass distributions

The LIGO Scientific Collaboration (2021) gives recent estimations on the possible shape of mass distribution for the more massive of the merging BHs (m1) based on so far detected BH-BH mergers. The proposed fit is a power-law mass function p ( m 1 ) ∝ m 1 − α $ p(m_1) \propto m_1^{-\alpha} $, with a break around 40 M⊙, with the first exponent for the masses m1 ≲ 40 M⊙, α 1 = 1 . 58 − 0.86 + 0.82 $ \alpha_1 = 1.58^{+0.82}_{-0.86} $ and the second exponent for m1 ≳ 40 M⊙, α 2 = 5 . 6 − 2.5 + 4.1 $ \alpha_2 = 5.6^{+4.1}_{-2.5} $.

We present the distributions of more massive, less massive BH masses, and total mass of BH-BH systems which merged at redshifts z <  2 for two tested models with revised RLOF physics. Figures 6 and 7 shows results for models M480.B and M481.B, respectively. We approximate more BH mass distributions (Figs. 6 and 7; top panels) with power-law function for an easy visual comparison with the LIGO/Virgo estimate.

thumbnail Fig. 6.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2 and O3. Results for model M480.B.

thumbnail Fig. 7.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2, and O3. Results for model M481.B.

The shape of more massive BH distribution for our model M480.B is in very good agreement with the shape of the BH-BH detections (Fig. 6, top panel). The distribution for this model may be divided as well into two parts, with the break around 45 M⊙. The first part, up to the break is characterized by a slow decrease corresponding to the exponent α1 ≈ 1.5 while after the break, the curve become much steeper with α2 ≈ 5.3 and the cut at ∼55 M⊙ (PSN, begining of the second mass gap). We note that this model distribution is a better match to the LIGO/Virgo estimate than recently published StarTrack results (see Fig. 15 for model M30.B in Abbott et al. (2020b): single power-law with index α = 3.6). The LIGO/Virgo total BH-BH system masses are within range of our simulated population of BH-BH mergers for model M480.B. The only exception is GW190521: the most massive detected BH-BH merger (Abbott et al. 2020b), which has both BHs (85 M⊙ + 66 M⊙) within our adopted pair-instability mass gap: 55−135 M⊙ (Belczynski et al. 2020). However, even such massive BHs can be possibly produced by massive stars if uncertainties on nuclear reaction rates and mixing in stellar interiors are taken into account (Belczynski 2020).

In the case of model M481.B, the more massive BH mass distribution is not a good match to LIGO/Virgo data. Although two power-law exponents (α1 ≈ 1.0 and α2 ≈ 5.0) are within LIGO/Virgo estimates, the break appears at quite a low mass Mbreak ≈ 15−20 M⊙. We note also the total BH-BH mass distribution ends at lower mass (Mcutoff = 100 M⊙) as compared with model M480.B (Mcutoff = 110 M⊙).

We also plot mass distributions for the model M380.B which represents our standard RLOF physics (Fig. 8). Similarly to M481.B, the break point of more massive BH mass distribution is at rather low BH masses, Mbreak ≈ 15−20 M⊙. However, in this case for masses m1 <  Mbreak the distribution is increasing with BH mass (∝M+2.7), which is inconsistent with the LIGO/Virgo estimates. For masses m1 >  Mbreak we fit power-law exponent α2 = 3.3. The maximum total mass (m1 + m2) of the BH-BH merger for model M380.B is 110 M⊙.

thumbnail Fig. 8.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2 and O3. Results for model M380.B (standard).

Our results indicate that modifications in RLOF physics approach may drastically influence the BH mass distribution of BH-BH mergers. Some approaches (M480.B) allow to match nicely the empirical data for more massive BH mass distribution while for others distribution shape is rather off (M380.B and M481.B).

5. Formation scenarios

5.1. BH-BH formation scenarios

In Table 4, we provide evolutionary scenarios which lead to the formation of local (z ∼ 0) BH-BH mergers for tested models with revised CE development criteria (M480.B and M481.B) and standard StarTrack criteria (M380.B). For each model the dominant formation scenario, which constitutes the highest merger rate density fraction (marked with the bold text) is different. In Figs. 9 and 10 we present diagrams with the example for the dominant BH-BH formation scenario in model M480.B (without CE phase), together with alternative evolution of the same system (the same initial conditions), but in models M380.B and M481.B. The Fig. 9 is an example of a low-mass BH-BH merger formation, for which significant natal kick (received after second SN explosion) plays an important role. The second example (Fig. 10) is massive BH-BH merger formation in which two BHs form via direct collapse.

thumbnail Fig. 9.

Typical local (z ≈ 0) BH-BH mergers formation scenarios for model M480.B together with the alternative evolution of the same system in models M380.B and M481.B. For all models we begin with the same initial conditions, massive system of 55 M⊙ primary and 35 M⊙ secondary on the orbit of 350 R⊙ and with metallicity Z = 0.0085. In all scenarios, the primary initiates TTMT when it leaves the main sequence and begins to expand in thermal timescale during the HG phase. After 4.8 Myr the primary finishes evolution and explodes as type Ib/Ic SN leaving behind a 11 M⊙ BH remnant. Next, after the secondary leaves main sequence, due to more restricted condition for CE development, in the cases of model M480.B and 481.B system goes through TTMT instead of CE phase as in standard M380.B model. During TTMT the secondary loses a substantial percentage of its mass (over 60%), together with the system orbital angular momentum. Therefore, the orbit tightens (by a factor of ∼4). After TTMT, in model 480.B, the system remains a BH-He system on the orbit of 150 R⊙. Due to the high natal kick after the second SN explosion, the orbital eccentricity significantly increases to e = 0.995, which allows the BH-BH system to merge in Hubble time. In standard (M380.B) scenario systems goes through CE phase with HG donor and we assume system a merger (see Sect. 2). In model M481.B due to modified condition for switch between the TTMT and nuclear-timescale stable MT, and the “safety” condition (see Sect. 3.2) after TTMT system also enters the CE phase with the HG donor and merges, leaving behind a single BH. The “He” in the diagram stands for a stripped helium core.

thumbnail Fig. 10.

Typical local (z ≈ 0) BH-BH mergers formation scenarios for model M480.B, together with the alternative evolution of the same system in models M380.B and M481.B. For all models, we begin with the same initial conditions, massive system of 145 M⊙ primary and 83 M⊙ secondary on the orbit of 263 R⊙ and with metallicity Z = 0.0001. In all scenarios, the primary initiates TTMT when it leaves main sequence and begins to expand in thermal timescale during the HG phase. When TTMT finishes, in case of models M380.B and M480.B secondary looses envelope. A BH with stripped helium core companion on the orbit of 350 R⊙ is left. In model M481.B due to modified condition for switch between the TTMT and nuclear-timescale stable MT, and the “safety” condition (read Sect. 3.2) when TTMT ends, the system goes through the CE phase. As the donor is a HG star, the CE ends with the system merger leaving behind an unevolved single star. After 3.3 Myr, in models M380.B and M340.B, the primary finishes evolution and collapse to a 20 M⊙ BH remnant. Next, after the secondary leaves main sequence, due to more restricted conditions for CE development, in cases of model M480.B system goes through TTMT instead of CE phase as in standard M380.B model. At the onset of TTMT, the secondary (donor) is over four times more massive than a BH companion. During the TTMT, the donor star looses ∼60 % of its mass together with the system orbital angular momentum. This leads to the significant orbit tightening by a factor of ∼12. Efficient orbital angular momentum loss during TTMT allows for the subsequent formation of 20 M⊙ BH and 30 M⊙ BH systems to merge in Hubble time. In the scenario for M380.B, the highly unequal mass system enters the CE phase with the HG donor. The binary system does not survive CE phase leaving behind a single BH. The “He” in the diagram stands for stripped helium core.

For model M480.B (revised CE development criteria and standard switch from TTMT to nuclear-timescale stable MT), the dominant BH-BH mergers formation scenario is without any CE phase. Similar channels were reported in many other recent works (van den Heuvel et al. 2017; Andrews et al. 2020; Zevin et al. 2021; Bavera et al. 2020; Marchant et al. 2021), as a possible formation scenario of LIGO/Virgo sources via isolated binary evolution. The typical evolution of binary system in model M480.B is as follow: When the primary (initially more massive star) leaves its main sequence and expands, system goes through a period of TTMT. We note that in some cases the donor star enters TTMT while it is at its HG. Survival of a binary through this phase is facilitated by high TTMT rate from the donor star. This removes most of (or entire) donor’s H-rich envelope and does not allow for a significant increase in donor radius, which, in turn, does not allow for development of CE phase and possible merger of the donor with its companion. After RLOF, the primary becomes a naked helium star. When the primary star completes its evolution, the first BH is formed. In the meantime, the secondary evolves off the main sequence, expands, and initiates a second phase of TTMT. The second RLOF, as well as the first one, is often initiated by a HG stars expanding on the thermal scales. If the system components have highly unequal masses (the donor mass is by a factor of 3 or more greater than the BH mass), the substantial mass loss from the system carries away also lots of the orbital angular momentum (even 90% comparing with the amount at the second RLOF onset). This leads to significant orbit tightening. Such a mechanism was reported previously, for instance, by van den Heuvel et al. (2017), Marchant et al. (2021). After MT, the secondary becomes naked helium star, and next it forms a second, typically less massive BH.

In all tested models of BHs with masses MBH <  10−15 M⊙, they may receive a natal kick at the time of their formation (see Sect. 2). The lower the BH mass, the greater probability of a high natal kick. Therefore, in the BH-BH formation scenario the second-born BH (usually the less massive one) may receive a significant kick. High natal kick may impart significant eccentricity to newly formed BH-BH system what decreases the inspiral timescale Tins ∼ (1 − e2)7/2, (Peters 1964), allowing such BH-BH system to merge in Hubble time. For example BH-BH system of 25 M⊙ and 10 M⊙ on the wide orbit of ∼900 R⊙ may merge in ∼10 Gyr if the eccentricity is as high as e = 0.992. For a comparison, in the case of circular orbit (e = 0) initial separation of the same BH-BH system would have to be less than ∼30 R⊙, for eccentricity e = 0.7 less than ∼45 R⊙, and for e = 0.9 less than ∼110 R⊙. The first stable RLOF phase in this scenario (M480.B) would be replaced in our standard model (M380.B) with the CE phase for many binaries as revised CE development criteria is more restrictive. Obviously, CE may lead in some cases to binary component merger and formation of a massive single stars. We note however, that even in our standard scenario dominant formation channel, the first RLOF is stable TTMT from HG star (compare dominant path 3 for model M380.B with dominant path 1 for model M480.B in Table 4). The second stable RLOF phase would be replaced by CE in the standard model. The CE during second RLOF is the main phase which leads the orbit size decrease for final BH-BH system to merge within Hubble time (see path 3 in Table 4 for model M380.B).

For Model M481.B (revised CE development criteria and modified condition for switch from TTMT to nuclear timescale MT), the main BH-BH formation channel also differs from the one in the standard Startrack approach (M380.B). Most of the evolution resembles the scenario described for M480.B: first TTMT initiated by the primary, then the first BH formation with a possible natal kick and the second TTMT initiated by the secondary. However, in model M481.B, due to a modified condition for switch from the TTMT to nuclear-timescale stable MT, and the “safety” condition (see Sect. 3.2), when TTMT ends, the system goes through the CE phase. If the system survives CE, the separation is reduced and the secondary’s envelope is ejected, leaving behind a naked helium core. Soon thereafter, in the CE phase (in some cases also already during CE), the secondary’s core collapse and the second (typically less massive) BH is formed.

5.2. BH-NS formation scenarios

In Table 5, we provide formation scenarios for local (z ≈ 0) BH-NS mergers for tested models with revised CE development criteria (M480.B and M481.B) and standard Startrack criteria (M380.B). For each model, the dominant formation scenario which constitutes the highest merger rate density fraction is marked with the bold text. For two models M380.B and M480.B, the dominant formation scenario is the same, however, for model M480.B, it constitutes over 70%, while for M380.B about 50% of total BH-NS merger rates. The evolutionary diagrams with the dominant formation scenarios are shown in Fig. 11.

Table 5.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different BH-NS formation scenarios.

thumbnail Fig. 11.

Typical local (z ≈ 0) BH-NS mergers formation scenarios for three tested models M380, M480, and M481 (for models M380 and M480 the scenario is the same). For all models, we begin with the same initial conditions: a wide binary system with 22 M⊙ primary and 11 M⊙ secondary on the orbit of ∼2600 R⊙ with metallicity Z = 0.01. Most of the evolutionary phases are common for all three models. When the primary leaves main sequence it expands, first as an HG star and then as core helium burning giant. The system goes through stable RLOF phase with TTMT and nuclear-timescale MT episodes. After MT the donor looses its envelope and soon explodes as SN. A low-mass BH (3.5 M⊙) is formed with the accompanying high natal kick which causes the system eccentricity increase to e = 0.83. Next, the secondary leaves main sequence and begins to expands. At the same time the orbit gets circularized by the tides. System goes thorough the CE phase, the secondary looses its envelope. After CE separation decreases to only 3 R⊙. Close system of BH and helium star is formed. As after CE the secondary is a main sequence or HG naked helium star, it expands initiating a stable RLOF. The system goes through another TTMT, which in model M481 ends with the second CE phase and formation of a very close (0.11 R⊙) BH-NS system. In the case of models M380 and M480, once TTMT ends, the system remains a close BH-He star binary and soon the secondary explodes forming a NS.

Most of the evolutionary phases are common for all three tested models. The evolutionary scenario transpires as follows: the once primary, initially more massive star, leaves the main sequence, it expands and initiates the first stable RLOF. System goes through the stable MT, after which the donor (primary) star looses its envelope becoming a stripped helium star. Next, the primary forms a BH through SN explosion or via direct collapse (depends on the final star mass). When the secondary star leaves its main sequence it expands and initiates unstable MT. The CE development leads to the secondary’s envelope removal and significant contraction of the system separation. Close system of BH and naked helium star is formed. Since, following the CE phase, the secondary is a main sequence or HG helium star, it still expands initiating a stable RLOF. At the end of the TTMT, the evolutionary scenarios for different models begin to diverge. In the case of models M380.B and M480.B, after the TTMT system is close to the BH-He binary. Next, a SN explosion takes place, followed by NS formation. In scenario M481.B, due to a modified condition for the switch between the TTMT and nuclear-timescale stable MT, and the “safety” condition (read Sect. 3.2), the system ends stable RLOF with CE phase. The orbit is tightened once again. During the second CE phase, the NS is formed.

6. Conclusions and discussion

In this work, we implement the revised RLOF stability diagram for binary systems with initially massive donors (MZAMS >  18 M⊙) into StarTrack population synthesis code. The revised diagram is based on MT simulations performed by Pavlovskii et al. (2017) using the 1D hydrodynamical stellar code MESA. They found that RLOF is stable over a much wider parameter space than previously thought if the MT is allowed to proceed until outflows from the outer Langrangian points develop. For one tested model, we also modify the conditions for switching between TTMT and nuclear timescale MT. We used the most updated version of StarTrack code to calculate local merger rate density and mass ratio distribution of the synthetic population of DCOs mergers.

We compare the results of the revised RLOF treatment with the standard StarTrack method. We also check if the new results are consistent with the recent estimates of LIGO/Virgo collaboration, comparing local merger rate density, the fraction of unequal mass mergers, and the shape of BH-BH mergers mass distribution.

We present the results of three models; our standard Star Track input physics model with delayed SN engine (M380.B) and two models with revised RLOF treatment (M480.B and M481.B). Model M480.B differs from our standard M380.B model only by the revised CE development criteria, while M481.B additionally includes the modified condition for the switch from TTMT to nuclear timescale MT.

Various changes related to assumptions about RLOF physics may strongly affect (sometimes in a non-intuitive way) the formation of DCO mergers: impacting merger rates and physical properties of merging objects. For the three presented models, we obtained not only different local merger rate densities or BH mergers mass distributions, but also the dominant evolutionary scenarios for BH-BH formation changes.

In model M480.B (with revised CE development criteria), the most common scenario for BH-BH merger formation is without any CE phase. In model M481.B (with revised CE development criteria and modified switch from TTMT to nuclear timescale MT), BH-BH dominant formation scenario involves CE that is preceded by short TTMT phase.

Among the tested models, we did not find one that would be fully consistent with the merger rates for NS-NS, BH-NS, and BH-BH given by LIGO/Virgo and at the same time would produce the reported more massive BH mass distribution for BH-BH mergers. While the shape of the more massive BH mass distribution in model M480.B fits two exponent values well as well as the break point given (The LIGO Scientific Collaboration 2021), the total rates for BH-BH mergers are too high by a factor of ∼2−3. Additionally, the fraction of unequal BH-BH mergers is not high enough. For Model M481.B, the total local rate density for all types of DCO mergers are aptly consistent with the most recent LIGO/Virgo estimates, but the shape of the more massive BH distribution seems too steep. Obviously, we could try to modify RLOF physics further or even change other parts of input physics (e.g., increasing the BH natal kicks to lower BH-BH merger rates) to fit LIGO/Virgo data for some models. However, the aim of this paper is only to demonstrate how uncertain ingredients of population synthesis codes may dramatically alter predictions.

We checked that the fraction of accreted mass during non-conservative RLOF (non-degenerate accretors, see Sect. 2) does not significantly influence the population of BH-BH and BH-NS mergers in the tested model M380.B. Beside the standard Startrack value fa = 0.5 , we tested models with values fa = 0.8 and fa = 0.2 for accreted mass fraction (the fraction 1 − fa of transferred mass is ejected from the system). We found that merger rates for BH-BH and BH-NS systems change about 3% for fa = 0.8 and about 20% for fa = 0.2 comparing to the model with fa = 0.5. Also, the mass and mass ratio distribution for three tested values are similar, with m1 (more massive of merging components) slightly shifted to the lower mass values in the model with the lowest accretion fraction fa = 0.2. Different adopted values of fa influence, however, the local merger rate density for NS-NS systems. In both fa = 0.2 and fa = 0.8 models, the rate for NS-NS systems decreased by a factor of ∼2 in comparison to the standard fa = 0.5 model. The drop in NS-NS rates is a complex mix of multiple overlapping physical processes, for instance different amounts of lost mass (together with orbital angular momentum) leads to different types of system separation after the first RLOF, which further leads to a different evolutionary type of the donor during the second RLOF (CE phase with HG donor ends with a system merger).

In earlier works based on StarTrack code, two main variants of CE treatment were usually considered (referred to as submodel A and B). Submodel B (the one presented in this work, see Sect. 2) assumes that all HG donor systems merge during the CE phase. In contrast, submodel A allows such binary systems to survive if the energy budget calculations allow. We decided to present only the results for submodel B for two main reasons: First, we checked that typically HG donors which initiate CE in StarTrack are only partially expanded stars, shortly after their main sequence. Such stars likely do not have a clear core-envelope boundary and, therefore, it seems more physical to assume their merger (as we do for main sequence stars). Second, the merger rates for DCOs, especially BH-BH systems in submodel A (see Belczynski et al. 2020) are often much larger (up to 30 times) than estimated by LIGO/Virgo.

The central message of our rapid population synthesis study is that the choice of our RLOF treatment produces noticeably different populations of GW sources. Our results highlight the need for caution when makes inferences from stellar and binary evolution models, as applied to LIGO/Virgo results. Similar cautionary notes are also found through more detailed evolutionary calculations with hydrodynamical stellar codes (Klencki et al. 2021; Decin 2020). If we suppose a population synthesis study is performed with variations of only natal kicks and CE efficiency, and the match to LIGO/Virgo rates is achieved for some specific CE efficiency and natal kicks assumption, it cannot be expected that we have already constrained these uncertain components of the input physics.

Acknowledgments

We would like to thank anonymous referee for their useful comments. K.B. and A.O. acknowledge support from the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050). N.I. acknowledges support from CRC program and funding from NSERC Discovery under Grant No. NSERC RGPIN-2019-04277.

References

  1. Abbott, B., Abbott, R., Abbott, T., et al. 2019a, Phys. Rev. X, 9, [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, ApJ, 882, L24 [NASA ADS] [CrossRef] [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019c, Phys. Rev. X, 9, 031040 [Google Scholar]
  4. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, ApJ, 896, L44 [Google Scholar]
  5. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, Phys. Rev. Lett., 125, 101102 [CrossRef] [Google Scholar]
  6. Andrews, J. J., Cronin, J., Kalogera, V., Berry, C., & Zezas, A. 2020, ApJ, 914, L32 [Google Scholar]
  7. Antonini, F., & Perets, H. B. 2012, ApJ, 757, 27 [NASA ADS] [CrossRef] [Google Scholar]
  8. Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77 [NASA ADS] [CrossRef] [Google Scholar]
  9. Arca-Sedda, M., & Capuzzo-Dolcetta, R. 2019, MNRAS, 483, 152 [NASA ADS] [CrossRef] [Google Scholar]
  10. Arca-Sedda, M., Li, G., & Kocsis, B. 2021, A&A, 650, A189 [EDP Sciences] [Google Scholar]
  11. Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bae, Y.-B., Kim, C., & Lee, H. M. 2014, MNRAS, 440, 2714 [NASA ADS] [CrossRef] [Google Scholar]
  13. Banerjee, S. 2018, MNRAS, 473, 909 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Begelman, M. C. 1979, MNRAS, 187, 237 [NASA ADS] [Google Scholar]
  16. Belczynski, K. 2020, ApJ, 905, L15 [Google Scholar]
  17. Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 [NASA ADS] [CrossRef] [Google Scholar]
  18. Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504 [NASA ADS] [CrossRef] [Google Scholar]
  19. Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
  20. Belczynski, K., Dominik, M., Bulik, T., et al. 2010, ApJ, 715, L138 [NASA ADS] [CrossRef] [Google Scholar]
  21. Belczynski, K., Wiktorowicz, G., Fryer, C. L., Holz, D. E., & Kalogera, V. 2012, ApJ, 757, 91 [NASA ADS] [CrossRef] [Google Scholar]
  22. Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016a, Nature, 534, 512 [Google Scholar]
  23. Belczynski, K., Heger, A., Gladysz, W., et al. 2016b, A&A, 594, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104 [CrossRef] [EDP Sciences] [Google Scholar]
  25. Benacquista, M. J., & Downing, J. M. B. 2013, Liv. Rev. Rel., 16, 4 [NASA ADS] [CrossRef] [Google Scholar]
  26. Bethe, H. A., & Brown, G. E. 1998, ApJ, 506, 780 [NASA ADS] [CrossRef] [Google Scholar]
  27. Bond, J. R., & Carr, B. J. 1984, MNRAS, 207, 585 [CrossRef] [Google Scholar]
  28. Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, ArXiv e-prints [arXiv:2103.02608] [Google Scholar]
  29. Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017, ApJ, 836, L26 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cherepashchuk, A. M., Postnov, K. A., & Belinski, A. A. 2019, MNRAS, 485, 2638 [CrossRef] [Google Scholar]
  31. Chruślińska, M., Jeřábková, T., Nelemans, G., & Yan, Z. 2020, A&A, 636, A10 [CrossRef] [EDP Sciences] [Google Scholar]
  32. Clayton, M., Podsiadlowski, P., Ivanova, N., & Justham, S. 2017, MNRAS, 470, 1788 [NASA ADS] [CrossRef] [Google Scholar]
  33. de Kool, M. 1990, ApJ, 358, 189 [NASA ADS] [CrossRef] [Google Scholar]
  34. de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545 [Google Scholar]
  35. Decin, L. 2020, ArXiv e-prints [arXiv:2011.13472] [Google Scholar]
  36. Dewi, J. D. M., & Tauris, T. M. 2000, A&A, 360, 1043 [NASA ADS] [Google Scholar]
  37. Di Carlo, U. N., Giacobbo, N., Mapelli, M., et al. 2019, MNRAS, 487, 2947 [NASA ADS] [CrossRef] [Google Scholar]
  38. Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52 [NASA ADS] [CrossRef] [Google Scholar]
  39. Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2010, MNRAS, 407, 1946 [NASA ADS] [CrossRef] [Google Scholar]
  40. Drozda, P., Belczynski, K., O’Shaughnessy, R., Bulik, T., & Fryer, C. L. 2020, ArXiv e-prints [arXiv:2009.06655] [Google Scholar]
  41. du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941 [Google Scholar]
  42. Eggleton, P. P. 2002, ApJ, 575, 1037 [NASA ADS] [CrossRef] [Google Scholar]
  43. Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302 [NASA ADS] [CrossRef] [Google Scholar]
  44. Fishbach, M., & Holz, D. E. 2020, ApJ, 891, L27 [CrossRef] [Google Scholar]
  45. Fragione, G., & Kocsis, B. 2019, MNRAS, 486, 4781 [CrossRef] [Google Scholar]
  46. Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, MNRAS, 488, 47 [CrossRef] [Google Scholar]
  47. Fragos, T., Andrews, J. J., Ramirez-Ruiz, E., et al. 2019, ApJ, 883, L45 [NASA ADS] [CrossRef] [Google Scholar]
  48. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [NASA ADS] [CrossRef] [Google Scholar]
  49. Ge, H., Hjellming, M. S., Webbink, R. F., Chen, X., & Han, Z. 2010, ApJ, 717, 724 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2015, ApJ, 812, 40 [NASA ADS] [CrossRef] [Google Scholar]
  51. Ge, H., Webbink, R. F., & Han, Z. 2020a, ApJS, 249, 9 [CrossRef] [Google Scholar]
  52. Ge, H., Webbink, R. F., Chen, X., & Han, Z. 2020b, ApJ, 899, 132 [CrossRef] [Google Scholar]
  53. Grichener, A., & Soker, N. 2021, ArXiv e-prints [arXiv:2101.05118] [Google Scholar]
  54. Gültekin, K., Miller, M. C., & Hamilton, D. P. 2004, ApJ, 616, 221 [NASA ADS] [CrossRef] [Google Scholar]
  55. Gültekin, K., Miller, M. C., & Hamilton, D. P. 2006, ApJ, 640, 156 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hainich, R., Oskinova, L. M., Shenar, T., et al. 2018, A&A, 609, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Hamers, A. S., Bar-Or, B., Petrovich, C., & Antonini, F. 2018, ApJ, 865, 2 [CrossRef] [Google Scholar]
  58. Han, Z., Podsiadlowski, P., & Eggleton, P. P. 1994, MNRAS, 270, 121 [NASA ADS] [CrossRef] [Google Scholar]
  59. Han, Z.-W., Ge, H.-W., Chen, X.-F., & Chen, H.-L. 2020, Res. Astron. Astrophys., 20, 161 [CrossRef] [Google Scholar]
  60. Hartwig, T., Volonteri, M., Bromm, V., et al. 2016, MNRAS, 460, L74 [NASA ADS] [CrossRef] [Google Scholar]
  61. Hoang, B.-M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, ApJ, 856, 140 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974 [NASA ADS] [CrossRef] [Google Scholar]
  63. Horvath, J. E., Rocha, L. S., Bernardo, A. L. C., de Avellar, M. G. B., & Valentim, R. 2020, ArXiv e-prints [arXiv:2011.08157] [Google Scholar]
  64. Howitt, G., Stevenson, S., Vigna-Gómez, A., et al. 2020, MNRAS, 492, 3229 [Google Scholar]
  65. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
  67. Hurley, J. R., Sippel, A. C., Tout, C. A., & Aarseth, S. J. 2016, MNRAS, 33 [Google Scholar]
  68. Iaconi, R., & De Marco, O. 2019, MNRAS, 490, 2550 [Google Scholar]
  69. Iaconi, R., Reichardt, T., Staff, J., et al. 2017, MNRAS, 464, 4028 [NASA ADS] [CrossRef] [Google Scholar]
  70. Ivanova, N., & Taam, R. E. 2004, ApJ, 601, 1058 [NASA ADS] [CrossRef] [Google Scholar]
  71. Ivanova, N., Justham, S., & Ricker, P. 2020, Common Envelope Evolution (IOP Publishing), 2514 [Google Scholar]
  72. Jones, D. 2020, Rev. Frontiers Mod. Astrophys., 123 [CrossRef] [Google Scholar]
  73. Kalogera, V., & Webbink, R. F. 1996, ApJ, 458, 301 [NASA ADS] [CrossRef] [Google Scholar]
  74. King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJ, 552, L109 [NASA ADS] [CrossRef] [Google Scholar]
  75. Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963 [NASA ADS] [CrossRef] [Google Scholar]
  76. Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54 [EDP Sciences] [Google Scholar]
  77. Kremer, K., Ye, C. S., Rui, N. Z., et al. 2020, ApJS, 247, 48 [CrossRef] [Google Scholar]
  78. Kroupa, P. 2002, Science, 295, 82 [Google Scholar]
  79. Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545 [Google Scholar]
  80. Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908 [Google Scholar]
  81. Lai, D., Rasio, F. A., & Shapiro, S. L. 1993, ApJS, 88, 205 [NASA ADS] [CrossRef] [Google Scholar]
  82. Leiner, E. M., & Geller, A. 2021, ApJ, 908, 229 [Google Scholar]
  83. Lipunov, V. M., Postnov, K. A., & Prokhorov, M. E. 1997, Astron. Lett., 23, 492 [NASA ADS] [Google Scholar]
  84. Liu, B., & Lai, D. 2018, ApJ, 863, 68 [CrossRef] [Google Scholar]
  85. Livio, M., & Soker, N. 1988, ApJ, 329, 764 [NASA ADS] [CrossRef] [Google Scholar]
  86. MacLeod, M., & Ramirez-Ruiz, E. 2015, ApJ, 803, 41 [NASA ADS] [CrossRef] [Google Scholar]
  87. MacLeod, M., Antoni, A., Murguia-Berthier, A., Macias, P., & Ramirez-Ruiz, E. 2017, ApJ, 838, 56 [NASA ADS] [CrossRef] [Google Scholar]
  88. Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
  89. Mandel, I., Müller, B., Riley, J., et al. 2021, MNRAS, 500, 1380 [CrossRef] [Google Scholar]
  90. Mapelli, M. 2016, MNRAS, 459, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  91. Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36 [Google Scholar]
  93. Marchant, P., Pappas, K. M. W., Gallegos-Garcia, M., et al. 2021, A&A, 650, A107 [CrossRef] [EDP Sciences] [Google Scholar]
  94. Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Meurs, E. J. A., & van den Heuvel, E. P. J. 1989, A&A, 226, 88 [NASA ADS] [Google Scholar]
  96. Miller, M. C., & Hamilton, D. P. 2002a, MNRAS, 330, 232 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  97. Miller, M. C., & Hamilton, D. P. 2002b, ApJ, 576, 894 [NASA ADS] [CrossRef] [Google Scholar]
  98. Misra, D., Fragos, T., Tauris, T. M., Zapartas, E., & Aguilera-Dena, D. R. 2020, A&A, 642, A174 [CrossRef] [EDP Sciences] [Google Scholar]
  99. Mondal, S., Belczyński, K., Wiktorowicz, G., Lasota, J.-P., & King, A. R. 2020, MNRAS, 491, 2747 [Google Scholar]
  100. Morawski, J., Giersz, M., Askar, A., & Belczynski, K. 2018, MNRAS, 481, 2168 [CrossRef] [Google Scholar]
  101. Nandez, J. L. A., & Ivanova, N. 2016, MNRAS, 460, 3992 [Google Scholar]
  102. Nandez, J. L. A., Ivanova, N., & Lombardi, J. C. J. 2015, MNRAS, 450, L39 [NASA ADS] [CrossRef] [Google Scholar]
  103. O’Leary, R. M., O’Shaughnessy, R., & Rasio, F. A. 2007, Phys. Rev. D, 76, 061504 [NASA ADS] [CrossRef] [Google Scholar]
  104. Olejak, A., Fishbach, M., Belczynski, K., et al. 2020, ApJ, 901, L39 [CrossRef] [Google Scholar]
  105. Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, eds. P. Eggleton, S. Mitton, & J. Whelan, IAU Symp., 73, 75 [Google Scholar]
  106. Paczyński, B., & Ziółkowski, J. 1968, Acta Astron., 18, 255 [Google Scholar]
  107. Passy, J.-C., De Marco, O., Fryer, C. L., et al. 2012, ApJ, 744, 52 [NASA ADS] [CrossRef] [Google Scholar]
  108. Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092 [NASA ADS] [CrossRef] [Google Scholar]
  109. Perna, R., Wang, Y.-H., Farr, W. M., Leigh, N., & Cantiello, M. 2019, ApJ, 878, L1 [NASA ADS] [CrossRef] [Google Scholar]
  110. Peters, P. C. 1964, PhD Thesis, California Institute of Technology [Google Scholar]
  111. Podsiadlowski, P., Ivanova, N., Justham, S., & Rappaport, S. 2010, MNRAS, 406, 840 [NASA ADS] [Google Scholar]
  112. Portegies Zwart, S. F., & Verbunt, F. 1996, A&A, 309, 179 [NASA ADS] [Google Scholar]
  113. Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  114. Ricker, P. M., & Taam, R. E. 2008, ApJ, 672, L41 [NASA ADS] [CrossRef] [Google Scholar]
  115. Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016, ApJ, 824, L8 [NASA ADS] [CrossRef] [Google Scholar]
  116. Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., et al. 2018, Phys. Rev. D, 98, 123005a [NASA ADS] [CrossRef] [Google Scholar]
  117. Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, ApJ, 676, 1162 [NASA ADS] [CrossRef] [Google Scholar]
  118. Sand, C., Ohlmann, S. T., Schneider, F. R. N., Pakmor, R., & Röpke, F. K. 2020, A&A, 644, A60 [CrossRef] [EDP Sciences] [Google Scholar]
  119. Santoliquido, F., Mapelli, M., Bouffanais, Y., et al. 2020, ApJ, 898, 152 [Google Scholar]
  120. Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021, MNRAS, 502, 4877 [CrossRef] [Google Scholar]
  121. Shiber, S., Iaconi, R., De Marco, O., & Soker, N. 2019, MNRAS, 488, 5615 [CrossRef] [Google Scholar]
  122. Silsbee, K., & Tremaine, S. 2017, ApJ, 836, 39 [NASA ADS] [CrossRef] [Google Scholar]
  123. Spera, M., Giacobbo, N., & Mapelli, M. 2016, Mem. Soc. Astron. It., 87, 575 [Google Scholar]
  124. Spera, M., Mapelli, M., Giacobbo, N., et al. 2019, MNRAS, 485, 889 [NASA ADS] [CrossRef] [Google Scholar]
  125. Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017, Nat. Commun., 8, 14906 [NASA ADS] [CrossRef] [Google Scholar]
  126. Tanikawa, A., Susa, H., Yoshida, T., Trani, A. A., & Kinugawa, T. 2021, ApJ, 910, 30 [CrossRef] [Google Scholar]
  127. The LIGO Scientific Collaboration& the Virgo Collaboration, 2020, ArXiv e-prints [arXiv:2004.08342] [Google Scholar]
  128. The LIGO Scientific Collaboration, the Virgo Collaboration (Abbott, R., et al.) 2021, ApJ, 913, L7 [CrossRef] [Google Scholar]
  129. Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675 [NASA ADS] [CrossRef] [Google Scholar]
  130. van den Heuvel, E. P. J. 1976, in Structure and Evolution of Close Binary Systems, eds. P. Eggleton, S. Mitton, & J. Whelan, IAU Symp., 73, 35 [NASA ADS] [CrossRef] [Google Scholar]
  131. van den Heuvel, E. P. J., Portegies Zwart, S. F., & de Mink, S. E. 2017, MNRAS, 471, 4256 [NASA ADS] [CrossRef] [Google Scholar]
  132. Vanbeveren, D. 1991, A&A, 252, 159 [NASA ADS] [Google Scholar]
  133. Vanbeveren, D., De Donder, E., Van Bever, J., Van Rensbergen, W., & De Loore, C. 1998, New Astron., 3, 443 [NASA ADS] [CrossRef] [Google Scholar]
  134. VanLandingham, J. H., Miller, M. C., Hamilton, D. P., & Richardson, D. C. 2016, ApJ, 828, 77 [NASA ADS] [CrossRef] [Google Scholar]
  135. Vigna-Gómez, A., MacLeod, M., Neijssel, C. J., et al. 2020, PASA, 37 [Google Scholar]
  136. Vinciguerra, S., Neijssel, C. J., Vigna-Gómez, A., et al. 2020, MNRAS, 498, 4705 [CrossRef] [Google Scholar]
  137. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  138. Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  139. Webbink, R. F. 1984, ApJ, 277, 355 [Google Scholar]
  140. Woosley, S. E. 2016, ApJ, 824, L10 [NASA ADS] [CrossRef] [Google Scholar]
  141. Woosley, S. E. 2017, ApJ, 836, 244 [Google Scholar]
  142. Xu, X.-J., & Li, X.-D. 2010, ApJ, 716, 114 [NASA ADS] [CrossRef] [Google Scholar]
  143. Zevin, M., Samsing, J., Rodriguez, C., Haster, C.-J., & Ramirez-Ruiz, E. 2019, ApJ, 871, 91 [NASA ADS] [CrossRef] [Google Scholar]
  144. Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152 [Google Scholar]
  145. Zuo, Z.-Y., & Li, X.-D. 2014, MNRAS, 442, 1980 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

List of tested evolutionary models.

In the text
Table 2.

Boundary radii values between the stable and unstable RLOF, values based on data from Pavlovskii et al. (2017).

In the text
Table 3.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different types of DCO systems.

In the text
Table 4.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different BH-BH formation scenarios.

In the text
Table 5.

Local (z ∼ 0) merger rate densities [Gpc−3 yr−1] for different BH-NS formation scenarios.

In the text

All Figures

thumbnail Fig. 1.

Hertzsprung-Russell diagram for massive stars with ZAMS masses: 20 M⊙, 32 M⊙, 50 M⊙, 80 M⊙, and 100 M⊙ at metallicity Z = 0.007; evolution of single star since main sequence to core helium burning. In the top panel, we marked stars cooler than log(Teff/K) < 3.73 with the red line and with the blue line, those stars with log(Teff/K) > 3.73. It is expected that the red line represents stars having convective envelopes, while the blue line is for radiative envelopes (Ivanova & Taam 2004; Belczynski et al. 2008). In the bottom panel, we marked with the blue line, the star main sequence and Hertzsprung gap phases and with the red line, the core helium burning phase. Evolutionary stages are defined as in formulas by Hurley et al. (2002) and Belczynski et al. (2008). With the grey lines, we marked the borders corresponding to stellar radii: 100 R⊙, 500 R⊙, and 1000 R⊙.
In the text
thumbnail Fig. 2.

Revised CE development criteria for stars with 10% of solar metallicity. On horizontal axis is the mass of the donor star. On the vertical axis is the radius of the donor. The circle and square points mark values of RS and RU, respectively (see Table 2). We fit lines to these data (see Eqs. (4)–(6)) to show regions of stable RLOF and unstable RLOF (CE). The color coded bar on the right denotes companion mass.
In the text
thumbnail Fig. 3.

Revised CE development criteria for stars with 100% of solar metallicity. Notation is the same as on Fig. 2.
In the text
thumbnail Fig. 4.

Distribution of mass ratio in BH-BH and BH-NS mergers (z ∼ 0). Top panel: BH-BH mergers, on the bottom panel BH-NS mergers. Blue line – results for standard CE development criteria with delayed SN engine (M380.B). Red line – results for revised CE development criteria (M480.B).
In the text
thumbnail Fig. 5.

Distribution of mass ratio in BH-BH and BH-NS mergers (z ∼ 0). On the top panel BH-BH mergers, on the bottom panel BH-NS mergers. Blue line – results for standard CE development criteria with delayed SN engine (M380.B). Red line – results for revised CE development criteria and revised stable RLOF treatment (M481.B).
In the text
thumbnail Fig. 6.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2 and O3. Results for model M480.B.
In the text
thumbnail Fig. 7.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2, and O3. Results for model M481.B.
In the text
thumbnail Fig. 8.

Intrinsic BH-BH merger (z <  2) mass distributions: More massive BH mass (top panel), less massive BH mass (middle panel), and total mass (bottom panel). We also show total mass estimates for LIGO/Virgo BH-BH mergers from O1, O2 and O3. Results for model M380.B (standard).
In the text
thumbnail Fig. 9.

Typical local (z ≈ 0) BH-BH mergers formation scenarios for model M480.B together with the alternative evolution of the same system in models M380.B and M481.B. For all models we begin with the same initial conditions, massive system of 55 M⊙ primary and 35 M⊙ secondary on the orbit of 350 R⊙ and with metallicity Z = 0.0085. In all scenarios, the primary initiates TTMT when it leaves the main sequence and begins to expand in thermal timescale during the HG phase. After 4.8 Myr the primary finishes evolution and explodes as type Ib/Ic SN leaving behind a 11 M⊙ BH remnant. Next, after the secondary leaves main sequence, due to more restricted condition for CE development, in the cases of model M480.B and 481.B system goes through TTMT instead of CE phase as in standard M380.B model. During TTMT the secondary loses a substantial percentage of its mass (over 60%), together with the system orbital angular momentum. Therefore, the orbit tightens (by a factor of ∼4). After TTMT, in model 480.B, the system remains a BH-He system on the orbit of 150 R⊙. Due to the high natal kick after the second SN explosion, the orbital eccentricity significantly increases to e = 0.995, which allows the BH-BH system to merge in Hubble time. In standard (M380.B) scenario systems goes through CE phase with HG donor and we assume system a merger (see Sect. 2). In model M481.B due to modified condition for switch between the TTMT and nuclear-timescale stable MT, and the “safety” condition (see Sect. 3.2) after TTMT system also enters the CE phase with the HG donor and merges, leaving behind a single BH. The “He” in the diagram stands for a stripped helium core.
In the text
thumbnail Fig. 10.

Typical local (z ≈ 0) BH-BH mergers formation scenarios for model M480.B, together with the alternative evolution of the same system in models M380.B and M481.B. For all models, we begin with the same initial conditions, massive system of 145 M⊙ primary and 83 M⊙ secondary on the orbit of 263 R⊙ and with metallicity Z = 0.0001. In all scenarios, the primary initiates TTMT when it leaves main sequence and begins to expand in thermal timescale during the HG phase. When TTMT finishes, in case of models M380.B and M480.B secondary looses envelope. A BH with stripped helium core companion on the orbit of 350 R⊙ is left. In model M481.B due to modified condition for switch between the TTMT and nuclear-timescale stable MT, and the “safety” condition (read Sect. 3.2) when TTMT ends, the system goes through the CE phase. As the donor is a HG star, the CE ends with the system merger leaving behind an unevolved single star. After 3.3 Myr, in models M380.B and M340.B, the primary finishes evolution and collapse to a 20 M⊙ BH remnant. Next, after the secondary leaves main sequence, due to more restricted conditions for CE development, in cases of model M480.B system goes through TTMT instead of CE phase as in standard M380.B model. At the onset of TTMT, the secondary (donor) is over four times more massive than a BH companion. During the TTMT, the donor star looses ∼60 % of its mass together with the system orbital angular momentum. This leads to the significant orbit tightening by a factor of ∼12. Efficient orbital angular momentum loss during TTMT allows for the subsequent formation of 20 M⊙ BH and 30 M⊙ BH systems to merge in Hubble time. In the scenario for M380.B, the highly unequal mass system enters the CE phase with the HG donor. The binary system does not survive CE phase leaving behind a single BH. The “He” in the diagram stands for stripped helium core.
In the text
thumbnail Fig. 11.

Typical local (z ≈ 0) BH-NS mergers formation scenarios for three tested models M380, M480, and M481 (for models M380 and M480 the scenario is the same). For all models, we begin with the same initial conditions: a wide binary system with 22 M⊙ primary and 11 M⊙ secondary on the orbit of ∼2600 R⊙ with metallicity Z = 0.01. Most of the evolutionary phases are common for all three models. When the primary leaves main sequence it expands, first as an HG star and then as core helium burning giant. The system goes through stable RLOF phase with TTMT and nuclear-timescale MT episodes. After MT the donor looses its envelope and soon explodes as SN. A low-mass BH (3.5 M⊙) is formed with the accompanying high natal kick which causes the system eccentricity increase to e = 0.83. Next, the secondary leaves main sequence and begins to expands. At the same time the orbit gets circularized by the tides. System goes thorough the CE phase, the secondary looses its envelope. After CE separation decreases to only 3 R⊙. Close system of BH and helium star is formed. As after CE the secondary is a main sequence or HG naked helium star, it expands initiating a stable RLOF. The system goes through another TTMT, which in model M481 ends with the second CE phase and formation of a very close (0.11 R⊙) BH-NS system. In the case of models M380 and M480, once TTMT ends, the system remains a close BH-He star binary and soon the secondary explodes forming a NS.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.