1. Introduction
A dynamically unstable mass transfer episode initiated by a post-main sequence (MS) donor is likely to lead to a common envelope episode (CEE), in which one star engulfs its companion and the binary spiral closer under the influence of drag forces (Paczynski Reference Paczynski, Eggleton, Mitton and Whelan1976). CEEs are proposed as a solution to the problem of how initially wide binaries, whose component stars may expand by tens to thousands of solar radii during their lifetime, become close binaries at later stages of evolution (van den Heuvel Reference van den Heuvel, Eggleton, Mitton and Whelan1976). Most evolutionary pathways leading to close compact binaries are expected to involve at least one CEE (Ivanova et al. Reference Ivanova2013a).
While CEEs are frequently invoked as a fundamental part of binary evolution, the detailed physics remain poorly understood (Iben & Livio Reference Iben and Livio1993; Ivanova et al. Reference Ivanova, Justham, Avendano Nandez and Lombardi2013b; Paczynski Reference Paczynski, Eggleton, Mitton and Whelan1976). There have been efforts in modelling and understanding the phase through hydrodynamic simulations, using Eulerian adaptive mesh refinement (Chamandy et al. Reference Chamandy2018; De et al. Reference De, MacLeod, Everson, Antoni, Mandel and Ramirez-Ruiz2019; Iaconi et al. Reference Iaconi2017, Reference Iaconi, De Marco, Passy and Staff2018; Li et al. Reference Li, Chang, Levin, Matzner and Armitage2020; L´opez-C´amara et al. Reference L´opez-C´amara, De Colle and Moreno M´endez2019; MacLeod & Ramirez-Ruiz Reference MacLeod and Ramirez-Ruiz2015; MacLeod et al. Reference MacLeod, Antoni, Murguia-Berthier, Macias and Ramirez-Ruiz2017b; Passy et al. Reference Passy2012; Ricker & Taam Reference Ricker and Taam2008, Reference Ricker and Taam2012; Sandquist et al. Reference Sandquist, Taam, Chen, Bodenheimer and Burkert1998; Shiber et al. 2019; Staff et al. Reference Staff, De Marco, Macdonald, Galaviz, Passy, Iaconi and Mac Low2016), moving meshes (Ohlmann et al. Reference Ohlmann, Röpke, Pakmor and Springel2016, Reference Ohlmann, Röpke, Pakmor and Springel2017; Prust & Chang Reference Prust and Chang2019), smoothed particle (Ivanova & Nandez Reference Ivanova and Nandez2016; Lombardi et al. Reference Lombardi, Proulx, Dooley, Theriault, Ivanova and Rasio2006; Nandez, Ivanova, & Lombardi Reference Nandez, Ivanova and Lombardi2015; Nandez & Ivanova Reference Nandez and Ivanova2016; Passy et al. Reference Passy2012; Rasio & Livio Reference Rasio and Livio1996; Reichardt et al. Reference Reichardt, De Marco, Iaconi, Tout and Price2019), particle-in-cell (Livio & Soker Reference Livio and Soker1988), and general relativistic (Cruz-Osorio & Rezzolla Reference Cruz-Osorio and Rezzolla2020) methods. Other approaches pursue detailed stellar modelling (Clayton et al. Reference Clayton, Podsiadlowski, Ivanova and Justham2017; Dewi & Tauris Reference Dewi and Tauris2000; Fragos et al. Reference Fragos, Andrews, Ramirez-Ruiz, Meynet, Kalogera, Taam and Zezas2019; Klencki et al. Reference Klencki, Nelemans, Istrate and Chruslinska2020; Kruckow et al. Reference Kruckow, Tauris, Langer, Sz´ecsi, Marchant and Podsiadlowski2016) or binary population synthesis (e.g., Andrews et al. Reference Andrews, Farr, Kalogera and Willems2015; Dewi, Podsiadlowski, & Sena Reference Dewi, Podsiadlowski and Sena2006; Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018; Nelemans et al. Reference Nelemans, Verbunt, Yungelson and Portegies Zwart2000; Tauris & Bailes Reference Tauris and Bailes1996; Vigna-Gómez et al. Reference Vigna-Gómez2018). There is currently no consensus on a thorough understanding of CEEs on all the relevant spatial and timescales.
Recent rapid population synthesis studies of double neutron star (DNS) populations have been partially motivated by the development of gravitational wave (GW) astronomy. Software tools such as StarTrack (Belczynski et al. Reference Belczynski2018; Chruslinska et al. Reference Chruslinska, Belczynski, Bulik and Gladysz2017, Reference Chruslinska, Belczynski, Klencki and Benacquista2018; Dominik et al. Reference Dominik2012), MOBSE (Giacobbo & Mapelli Reference Giacobbo and Mapelli2018, Reference Giacobbo and Mapelli2019a,b, Reference Giacobbo and Mapelli2020), ComBinE (Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018), and Compact Object Mergers: Population Astrophysics and Statistics (COMPAS) (Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020; Vigna-Gómez et al. Reference Vigna-Gómez2018) have been used to explore synthetic DNS populations in detail. Most of those studies focus on predicting or matching the observed DNS merger rate, either by investigating different parameterisations of the physics or varying the parameters within the models. In particular, all of the aforementioned population synthesis codes follow a similar simplified treatment of the common envelope (CE) phase.
In this paper, we present a detailed analysis of the CEEs that most merging DNSs are believed to experience at some point during their formation (Andrews et al. Reference Andrews, Farr, Kalogera and Willems2015; Belczynski, Kalogera, & Bulik Reference Belczynski, Kalogera and Bulik2002; Belczynski et al. Reference Belczynski2018; Bhattacharya & van den Heuvel Reference Bhattacharya and van den Heuvel1991; Dewi & Pols Reference Dewi and Pols2003; Dewi, Podsiadlowski, & Pols Reference Dewi, Podsiadlowski and Pols2005; Ivanova et al. Reference Ivanova, Belczynski, Kalogera, Rasio and Taam2003; Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018; Tauris & van den Heuvel Reference Tauris and van den Heuvel2006; Tauris et al. Reference Tauris2017; Vigna-Gómez et al. Reference Vigna-Gómez2018). Dominik et al. (Reference Dominik2012) previously used rapid population synthesis to study the relationship between CEE and DNS merger rates. In this study, we focus our attention on the properties at the onset of the Roche-lobe overflow (RLOF) episode leading to the CEE. We consider both long-period non-merging DNSs as well as short-period merging DNSs. We propose these distributions of binary properties as initial conditions for detailed studies of CEEs. We provide the results of this study in the form of a publicly available catalogueFootnote a.
We examine the properties of binaries unaffected by external dynamical interactions that experience CEEs on their way to forming DNS systems. We briefly discuss CEE leading to DNS formation in the context of generating some of the brightest luminous red novae, which may be signatures of CE ejections (Blagorodnova et al. Reference Blagorodnova2017; Howitt et al. Reference Howitt2020; Ivanova et al. Reference Ivanova, Justham, Avendano Nandez and Lombardi2013b; Pastorello et al. Reference Pastorello2019; MacLeod et al. Reference MacLeod, Macias, Ramirez-Ruiz, Grindlay, Batta and Montes2017a), and Be X-ray binaries (Vinciguerra et al. Reference Vinciguerra2020).
This paper is structured in the following way. Section 2 describes the initial distributions and relevant physical parameterisations used in rapid population synthesis. Section 3 presents the results of our study, particularly Hertzsprung–Russell (HR) diagrams displaying different properties of the systems, as well as their distributions. Section 4 discusses the results and some of the caveats. Finally, Section 5 summarises and presents the conclusions of this work.
2. Population synthesis model
We characterise CEEs with the rapid population synthesis element of the COMPAS suiteFootnote b (Barrett et al. Reference Barrett2018; Neijssel et al. Reference Neijssel2019; Stevenson et al. Reference Stevenson, Vigna-Gómez, Mandel, Barrett, Neijssel, Perkins and de Mink2017; Vigna-Gómez et al. Reference Vigna-Gómez2018). Rapid population synthesis relies on simplified methods and parameterisations in order to simulate a single binary from the zero-age main sequence (ZAMS) until stellar merger, binary disruption, or double compact object (DCO) formation. This approach relies on subsecond evolution of a single binary in order to generate a large population within hours using a single processor.
In COMPAS, an initial binary is defined as a gravitationally bound system completely specified by its metallicity, component masses, separation, and eccentricity at the ZAMS. We assume that our binaries have solar metallicity 
$Z=Z_{\odot}=0.0142$
(Asplund et al. Reference Asplund, Grevesse, Sauval and Scott2009). The mass of the primary (
$m_{\textrm{1}}$
), that is, the more massive star in the binary at birth, is drawn from the initial mass function 
$dN/dm_1 \propto m_{1}^{-2.3}$
(Kroupa Reference Kroupa2001) sampled between 
$5 \leq m_1/\textrm{M}_{\odot} \leq 100$
. The mass of the secondary (
$m_2$
) is obtained by drawing from a flat distribution in mass ratio (
$q_{\textrm{ZAMS}}=m_2/m_1$
) in the form 
$dN/dq \propto 1$
with 
$0.1 < q_{\textrm{ZAMS}} \leq 1$
(Sana et al. Reference Sana2012). The initial separation is drawn from a flat-in-the-log distribution in the form 
$dN/da \propto a^{-1}$
with 
$0.01 < a_{\textrm{ZAMS}}/\textrm{AU} < 1000$
(Öpik Reference Öpik1924). We assume that all our binaries have zero eccentricity at formation (the validity of this assumption is discussed in Section 4.4.3).
2.1. Adaptive importance sampling
COMPAS originally relied on Monte Carlo sampling from the birth distributions described above. However, this becomes computationally expensive when studying rare events.
In order to efficiently sample the parameter space leading to DNS formation, we adopt STROOPWAFEL as implemented in COMPAS (Broekgaarden et al. Reference Broekgaarden2019). STROOPWAFEL is an adaptive importance sampling (AIS) algorithm designed to improve the efficiency of sampling of unusual astrophysical events. The use of AIS increases the fraction of DNSs per number of binaries simulated by 
$\sim\! 2$
orders of magnitude with respect to regular Monte Carlo sampling. After sampling from a distribution designed to increase DNS yield, the binaries are reweighted by the ratio of the desired probability distribution of initial conditions to the actual sampling probability distribution.
We use bootstrapping to estimate the sampling uncertainty. We randomly resample each model population with replacement in order to generate a bootstrapped distribution. We perform this process 
$N=100$
times to get a 10% accuracy of the bootstrapped standard deviation. We calculate and report the standard deviation of the bootstrapped distributions as 
$1\sigma$
error bars.
2.2. Underlying physics
We mostly follow the physical model as presented in Vigna-Gómez et al. (Reference Vigna-Gómez2018). However, we highlight key aspects of the model that are particularly relevant for this work, along with a few non-trivial changes to the code. See also Appendix A for a summary and additional details on the set-up.
1. We approximate the Roche lobe radius following the fitting formula provided by Eggleton (Reference Eggleton1983) in the form:
(1)where
\begin{equation} \frac{R_{\textrm{RL}}}{a_{\textrm{p}}}= \frac{0.49q_{\textrm{RL}}^{2/3}}{0.6q_{\textrm{RL}}^{2/3}+ \textrm{ln}(1+q_{\textrm{RL}}^{1/3})}, 0<q_{\textrm{RL}}<\infty, \end{equation}
$R_{\textrm{RL}}$
is the effective Roche lobe radius of the donor,
$a_{\textrm{p}}=a(1-e)$
is the periastron, a and e are the semi-major axis and eccentricity, respectively,
$q_{\textrm{RL}}$
is the mass ratio;
$q_{\textrm{RL}}=m_{\textrm{donor}}/m_{\textrm{comp}}$
, with
$m_{\textrm{donor}}$
and
$m_{\textrm{comp}}$
being the mass of the donor and companion star, respectively. RLOF will occur once
$R_{\textrm{donor}} \geq R_{\textrm{RL}}$
, where
$R_{\textrm{donor}}$
is the radius of the donor.2. We use the properties of the system at the onset of RLOF in order to determine whether the mass transfer episode leads to a CEEs. Dynamical stability is determined by comparing the response of the radius of the donor to (adiabatic) mass loss to the response of the Roche lobe radius to mass transfer. This is done using the mass-radius exponent
(2)where the subscript ‘i’ represents either the mass-radius exponent for the donor (\begin{equation} \zeta_{\textrm{i}}=\frac{d\log R_{\textrm{i}}}{d\log m_{\textrm{donor}}}, \end{equation}
$\zeta_{\textrm{donor}}$
) or for the Roche lobe (
$\zeta_{\textrm{RL}}$
). We assume that
(3)leads to a CEE. Inspired by Ge et al. (Reference Ge, Webbink, Chen and Han2015), for MS donors, we assume\begin{equation} \zeta_{\textrm{donor}}
< \zeta_{\textrm{RL}} \end{equation}
$\zeta_{\textrm{donor}}=2.0$
; for Hertzsprung gap (HG) donors, we assume
$\zeta_{\textrm{donor}}=6.5$
. For post-helium ignition phases in which the donor still has a hydrogen envelope, we follow Soberman et al. (Reference Soberman, Phinney and van den Heuvel1997). All mass transfer episodes from stripped post-helium ignition stars, that is, case BB mass transfer (Delgado & Thomas Reference Delgado and Thomas1981; Dewi et al. Reference Dewi, Pols, Savonije and van den Heuvel2002; Dewi & Pols Reference Dewi and Pols2003) neutron star (NS) are assumed to be dynamically stable. For more details and discussion, see Tauris et al. (Reference Tauris, Langer and Podsiadlowski2015) and Vigna-Gómez et al. (Reference Vigna-Gómez2018).3. We deviate from Stevenson et al. (Reference Stevenson, Vigna-Gómez, Mandel, Barrett, Neijssel, Perkins and de Mink2017) and Vigna-Gómez et al. (Reference Vigna-Gómez2018) by allowing MS accretors to survive a CEE. Previously, any MS accretor was mistakenly assumed to imminently lead to a stellar merger. We now treat MS accretors just like any other stellar type. This does not have any effect on the COMPAS DNS population, as there are no dynamically unstable mass transfer phases with MS accretors leading to DNS formation (see discussion on formation history in Section 3.1).
4. We follow Kool (Reference de Kool1990) in the parameterisation of the binding energy (
$E_{\textrm{bind}}$
) of the donor star’s envelope (
$m_{\textrm{donor,env}}$
) given as:
(4)where G is the gravitational constant and\begin{equation} E_{\textrm{bind}}=\frac{-Gm_{\textrm{donor}}m_{\textrm{donor,env}}}{\lambda R_{\textrm{donor}}}, \end{equation}
$\lambda$
is a numerical factor that parameterises the binding energy.5. For the value of the
$\lambda$
parameter, we follow the fitting formulae from detailed stellar models as calculated by Xu & Li (Reference Xu and Li2010a,b). This
$\lambda$
, originally referred to as
$\lambda_{\textrm{b}}$
, includes internal energy and is implemented in the same way as
$\lambda_{\textrm{Nanjing}}$
in StarTrack (Dominik et al. Reference Dominik2012). Additionally, we fixed a bug which underestimated the binding energy of the envelope. We discuss the effect that has on the DNS population in Section 3.2.6. We use the
$\alpha\lambda$
-formalism (de Kool Reference de Kool1990; Webbink Reference Webbink1984) to determine the post-CEE orbit, with
$\alpha=1$
in all of our CEE.7. We use the Fryer et al. (Reference Fryer, Belczynski, Wiktorowicz, Dominik, Kalogera and Holz2012) delayed supernova (SN) remnant mass prescription, which was the preferred model from Vigna-Gómez et al. (Reference Vigna-Gómez2018). This prescription allows for a continuous (gravitational) remnant mass distribution between NSs and black holes (BHs), with a transition point at
$2.5\ \textrm{M}_{\odot}$
.8. The remnant mass of NSs with large baryonic mass previously only accounted for neutrino mass loss instead of an actual equation of state. This has now been corrected. This only affects NSs with remnant masses larger than
${\approx} 2.15\ \textrm{M}_{\odot}$
.
2.3. Tidal timescales
Mass transfer episodes occur in close binaries that experienced tidal interactions. The details of these tidal interactions are sensitive to the properties of the binary and the structure of the envelope of the tidally distorted stars, either radiative or convective.
The equilibrium tide refers to viscous dissipation in a star that is only weakly perturbed away from the shape that it would have in equilibrium (Zahn Reference Zahn1977). Meanwhile, the dynamical tide (Zahn Reference Zahn1975) refers to the excitation of multiple internal modes of a star in a time-varying gravitational potential; when these oscillatory modes are damped, orbital energy is lost to thermal energy (Eggleton, Kiseleva, & Hut Reference Eggleton, Kiseleva and Hut1998; Moe & Kratter Reference Moe and Kratter2018). Tidal evolution tends to align and synchronise the component spins with the direction of the orbital angular momentum vector and circularise the binary (Counselman 1973; Zahn Reference Zahn, Goupil, Zahn and Publications Series2008).
There are numerous uncertainties in tidal evolution. For example, the role of eccentricity is an active field of research. Heartbeat stars are eccentric binaries with close periastron passage which experience tidal excitation of different oscillatory modes (see Shporer et al. Reference Shporer2016 and references therein). Eccentric systems may also experience resonance locking, which occurs when a particular tidal harmonic resonates with a stellar oscillation mode; this enhances the efficiency of tidal dissipation (Witte & Savonije Reference Witte and Savonije1999a,b). The high-eccentricity regime, previously studied in the parabolic (and chaotic) limit (Mardling Reference Mardling1995a,b), has recently being revisited in the context of both dynamical (Vick & Lai Reference Vick and Lai2018) and equilibrium tides (Vick & Lai Reference Vick and Lai2020). There are uncertainties in the low-eccentricity regime, for example, there is a range of parameterisations for the equilibrium tide, such as the weak friction approximation, turbulent viscosity, and fast tides (Zahn Reference Zahn, Goupil, Zahn and Publications Series2008).
Here, we make several simplifying approximations for the synchronisation and circularisation timescales, 
$\tau_{\textrm{sync}}$
and 
$\tau_{\textrm{circ}}$
respectively, in order to parameterise the tidal evolution of the system.
We assume that the equilibrium tide operates on all stars with a convective envelope, regardless of the binary eccentricity. We use the equilibrium tide description in the weak friction model as described by Hut (Reference Hut1981) and implemented by Hurley, Tout, & Pols (Reference Hurley, Tout and Pols2002), although this may not be accurate for high-eccentricity systems (but see Vick & Lai Reference Vick and Lai2020). Since the equilibrium tide is generally a more efficient energy transport/dissipation mechanism than the dynamical tide for stars with convective envelopes, we ignore the contribution of the latter. Our equilibrium tide model is summarised in Section 2.3.1.
We apply the dynamical tide only to stars with a radiative envelope. In Section 2.3.2, we present our implementation of the dynamical tide following Zahn (Reference Zahn1977), as used in Hurley et al. (Reference Hurley, Tout and Pols2002).
2.3.1. The equilibrium tide for stars with convective envelopes
Under the equilibrium tide, the synchronisation and circularisation evolution equations for tides acting on a star of mass 
$m_{\textrm{tide}}$
from a companion star with mass 
$m_{\textrm{comp}}$
are
\begin{equation} \begin{split} \frac{d\Omega{\textrm{spin}}}{dt}= &\ 3\bigg( \frac{k}{\tau_{\textrm{tide}}} \bigg) \frac{q^2}{r_{\textrm{g}}^2} \bigg( \frac{R_{\textrm{tide}}}{a} \bigg)^6 \frac{\Omega_{\textrm{orb}}}{(1-e^2)^6} \\ &\times \bigg[ f_2(e^2)-(1-e^2)^{3/2}f_5(e^2)\frac{\Omega_{\textrm{spin}}}{\Omega_{\textrm{orb}}} \bigg], \end{split}\end{equation}and
\begin{equation} \begin{split} \frac{de}{dt}= &-27\bigg( \frac{k}{\tau_{\textrm{tide}}} \bigg) q(1+q) \bigg( \frac{R_{\textrm{tide}}}{a} \bigg)^8 \frac{e}{(1-e^2)^{13/2}} \\ &\times \bigg[ f_3(e^2)-\frac{11}{18}(1-e^2)^{3/2}f_4(e^2)\frac{\Omega_{\textrm{spin}}}{\Omega_{\textrm{orb}}} \bigg], \end{split}\end{equation}where 
$f_{n}(e^2)$
are polynomial expressions given by Hut (Reference Hut1981). The structure of the tidally deformed star is parameterised by k, which is the apsidal motion constant (Lecar, Wheeler, & McKee Reference Lecar, Wheeler and McKee1976) and the intrinsic tidal timescale (
$\tau_{\textrm{tide}}$
), usually associated with viscous dissipation (Zahn Reference Zahn1977). We follow Hurley et al. (Reference Hurley, Tout and Pols2002) in the calculation of the 
$(k/\tau_{\textrm{tide}})$
factor, which depends on the evolutionary stage and structure of the star. The mass ratio is defined as 
$q=m_{\textrm{comp}}/m_{\textrm{tide}}=1/q_{\textrm{RL}}$
and the gyration radius as 
$r_{\textrm{g}}=\sqrt{I_{\textrm{tide}}/(m_{\textrm{tide}}R_{\textrm{tide}}^2)}$
, where 
$I_{\textrm{tide}}$
and 
$R_{\textrm{tide}}$
are the moment of inertia and the radius of the tidally deformed star, respectively. The mean orbital velocity and the donor spin angular velocity are denoted by 
$\Omega_{\textrm{orb}}$
and 
$\Omega_{\textrm{spin}}$
, respectively.
Given that 
$a>R_{\textrm{tide}}$
, for a non-synchronous eccentric binary, we expect synchronisation to be faster than circularisation. If we assume that the system is synchronous (
$\Omega_{\textrm{orb}}=\Omega_{\textrm{spin}}$
), we simplify Equation (6) and estimate the circularisation timescale as
\begin{equation} \begin{split} \tau_{\textrm{circ}} = &\ -\frac{e}{de/dt} \\ = &\ \bigg\{ 27\bigg( \frac{k}{\tau_{\textrm{tide}}} \bigg) q(1+q) \bigg( \frac{R_{\textrm{tide}}}{a} \bigg)^8 \frac{1}{(1-e^2)^{13/2}} \\ &\ \times \bigg[ f_3(e^2)-\frac{11}{18}(1-e^2)^{3/2}f_4(e^2) \bigg] \bigg\}^{-1}. \end{split}\end{equation}2.3.2. The dynamical tide for stars with radiative envelopes
Following the derivation by Zahn (Reference Zahn1977), we can write the synchronisation and circularisation timescales for the dynamical tide as
\begin{equation} \begin{split} \tau_{\textrm{sync}} = &\ 52^{-5/3} \bigg( \frac{R_{\textrm{tide}}^3}{Gm_{\textrm{tide}}} \bigg)^{1/2} \frac{r_{\textrm{g}}^2}{q^2} (1+q)^{-5/6} \\ &\ \times E_{2}^{-1} \bigg( \frac{D}{R_{\textrm{tide}}} \bigg)^{17/2} \end{split}\end{equation}and
\begin{equation} \tau_{\textrm{circ}} = \frac{2}{21} \bigg( \frac{R_{\textrm{tide}}^3}{Gm_{\textrm{tide}}} \bigg)^{1/2} \frac{(1+q)^{-11/6}}{q} E_{2}^{-1} \bigg( \frac{D}{R_{\textrm{tide}}} \bigg)^{21/2},\end{equation}where 
$E_2=1.592\times 10^{-9} (M/\textrm{M}_{\odot})^{2.84}$
is a second-order tidal coefficient as fitted by Hurley et al. (Reference Hurley, Tout and Pols2002) from the values given by Zahn (Reference Zahn1975), under the assumption (violated for some of the systems we consider) that close binaries are nearly circular. For the dynamical tide, we set the tidal separation (D) to be the semilatus rectum 
$D=a(1-e^2)$
. This corresponds to the conservation of orbital angular momentum 
$J_{\textrm{orb}} \propto \sqrt{a(1-e^2)}$
. This assumption may lead us to underestimate the circularisation timescale for stars with radiative envelopes in highly eccentric orbits.
The dynamical tide is much less efficient than the equilibrium tide for virtually all binaries; therefore, we ignore the contribution of dynamical tides for convective envelope stars, even though they are active along with equilibrium tides.
Given the uncertainties in tidal circularisation efficiency, we do not include tides in dynamical binary evolution. Instead, we evolve binaries without the impact of tides, then estimate whether tides would have been able to circularise the binary prior to the onset of RLOF leading to a CEE as described below.
2.3.3. Radial expansion timescale
The strong dependence of the tidal timescales on 
$R_{\textrm{tide}}/a$
means that tides only become efficient when the star expands to within a factor of a few of the binary separation. Therefore, the rate of expansion of the star, which depends on the stage of stellar evolution, plays a key role in determining the efficiency of circularisation: the binary can circularise only if the circularisation timescale of an eccentric binary is shorter than the star’s radial expansion timescale. We define this radial expansion timescale as the radial e-folding time 
$\tau_{\textrm{radial}} \equiv dt/d\log R$
. This is computed by evaluating the local derivatives within the fitting formulae of Hurley, Pols & Tout (Reference Hurley, Pols and Tout2000) to the detailed stellar models from Pols et al. (Reference Pols, Schröder, Hurley, Tout and Eggleton1998).
2.3.4. Uncertainties in timescales
The timescales defined here, rather than fully accurate descriptions of tidal evolution, are used as order of magnitude estimates to analyse the overall properties of the population. Tidal timescales have significant uncertainties, including in the treatment of the dominant dissipation mechanism (e.g., weak friction approximation, turbulent convection, fast tides) and their parameterisation (Zahn Reference Zahn, Goupil, Zahn and Publications Series2008) and implementation (Hurley et al. Reference Hurley, Tout and Pols2002). Siess et al. (Reference Siess, Izzard, Davis and Deschamps2013) noted the problem with the 
$E_2$
fit being commonly misused, both via interpolation and extrapolation of stars above 
$20\ \textrm{M}_{\odot}$
(see also the alternative approach of Kushnir et al. Reference Kushnir, Zaldarriaga, Kollmeier and Waldman2017). The calculation of k and 
$\tau_{\textrm{tide}}$
follows Hurley et al. (Reference Hurley, Pols and Tout2000) and is uncertain for massive stars. For the radial expansion timescale, the fitting formulae we use are not accurate in representing the evolution of the star on thermal or dynamical timescales. These formulae also miss detailed information about the evolution of, for example, the size of the convective envelope. Additionally, they are not accurate in representing the effect of mass loss and mass gain.
3. Results
We present the results of the synthetic population of binaries which become DNSs. We focus our attention on the properties of the systems at the onset of the CEE. If a donor star experiences RLOF, leading to a dynamically unstable mass transfer episode, the system is classified as experiencing a CEE. In that case, we report the properties of the system at the moment of RLOF. We do not resolve the details of the CEE, such as the possible delayed onset of the dynamical inspiral phase. Given that we are interested in DNS progenitors, all of these CEEs will, by selection, experience a successful ejection of the envelope, that is, no stellar mergers are reported in this study. All the data presented in this work are available at https://zenodo.org/record/3593843 (Vigna-Gómez Reference Vigna-Gómez2019).
Our synthetic dataset contains about 1 000 000 binaries evolved using COMPAS. Out of all the simulated binaries, targeted at DNS-forming systems (see Section 2.1), there are 15 201 CEEs leading to DNS formation. These provide a far more accurate sampling of the 
${\approx}$
365 systems that would be expected for 86 000 000 
$\textrm{M}_{\odot}$
of star-forming mass sampled from the initial conditions. For simplicity, we assume 100% binarity a priori. Nevertheless, given our assumed separation distribution that is capped at 1000 AU, 10% of our systems never experience any mass transfer episode, resulting in two effectively single stars. While DNSs are believed to form in different environments, several studies have shown that metallicity does not play a large role in DNS properties, unlike binary BH or BH–NS formation (Dominik et al. Reference Dominik2012; Giacobbo & Mapelli Reference Giacobbo and Mapelli2018; Neijssel et al. Reference Neijssel2019; Vigna-Gómez et al. Reference Vigna-Gómez2018).
The results section is structured as follows. Section 3.1 discusses the two dominant formation channels in our model, that is, the evolutionary history of the binary from ZAMS to DNS formation. In Section 3.2, we present a comparison with the results from Vigna-Gómez et al. (Reference Vigna-Gómez2018). In Section 3.3, we describe the way main results are presented. In Section 3.4, we report the properties of the donor. In Section 3.5, we report the properties of the binary, in particular the orbital properties. Finally, in Section 3.6, we present and report the tidal circularisation timescales.
Schematic representation of DNS formation channels as described in Section 3.1. Top: Channel I is the dominant formation channel for DNS systems, as well as the most common formation channel in the literature (see, e.g., Tauris et al. Reference Tauris2017 and references therein). Bottom: formation Channel II distinguished by an early double-core CE phase. Acronyms as defined in text. Credit: T. Rebagliato.
3.1. Formation channels of DNS systems
Two common evolutionary pathways leading to the formation of DNS from isolated binary evolution are identified in the literature (Bhattacharya & van den Heuvel Reference Bhattacharya and van den Heuvel1991; Tauris & van den Heuvel Reference Tauris and van den Heuvel2006; Tauris et al. Reference Tauris2017). Following Vigna-Gómez et al. (Reference Vigna-Gómez2018), we refer to these formation channels as Channel I and Channel II.
Channel I is illustrated in the top panel of Figure 1 and proceeds in the following way:
1. A post-MS primary engages in stable mass transfer onto a MS secondary.
2. The primary, now stripped, continues its evolution as a naked helium star until it explodes in a SN, leaving a NS remnant in a bound orbit with a MS companion.
3. The secondary evolves off the MS, expanding and engaging in a CEE with the NS accretor.
4. After successfully ejecting the envelope, and hardening the orbit, the secondary becomes a naked helium star.
5. The stripped post-helium-burning secondary engages in highly non-conservative stable (case BB) mass transfer onto the NS companion.
6. After being stripped of its helium envelope, the ultra-stripped secondary (Tauris et al. Reference Tauris, Langer, Moriya, Podsiadlowski, Yoon and Blinnikov2013, Reference Tauris, Langer and Podsiadlowski2015) continues its evolution until it explodes as an ultra-stripped SN (USSN), forming a DNS.
In Channel I, the CEE may occur while the donor is crossing the HG, that is, between the end of the MS and the start of the core helium burning (CHeB) phase. Rapid population synthesis modelling of CEEs sometimes parameterise these donors in two possible outcomes: optimistic and pessimistic (Dominik et al. Reference Dominik2012). The optimistic approach assumes the donor has a clear core/envelope separation and that, as a result, the two stellar cores can potentially remove the CE, allowing the binary to survive the CEE. Throughout this paper, we assume the optimistic approach unless stated otherwise. The pessimistic approach assumes that dynamically unstable mass transfer from a HG donor leads imminently to a merger. The pessimistic approach results in 4% of potential DNS candidates merging before DCO formation.
Channel II is illustrated in the bottom panel of Figure 1 and proceeds in the following way:
1. A dynamically unstable mass transfer episode leads to a CEE when the primary and the secondary are both post-MS star. During this CEE, both stars have a clear core-envelope separation, and they engage in what is referred to in the literature as a double-core CEE (Brown Reference Brown1995; Dewi et al. Reference Dewi, Podsiadlowski and Sena2006; Justham, Podsiadlowski, & Han Reference Justham, Podsiadlowski and Han2011). For these binaries, evolutionary timescales are quite similar, with a minimum and mean mass ratio of
${\approx}0.93$
and
${\approx}0.97$
, respectively, consistent with high-mass and low-mass solar metallicity values reported in Dewi et al. (Reference Dewi, Podsiadlowski and Sena2006). During this double-core CEE, both stars are stripped and become naked helium stars.2. The stripped post-helium-burning primary engages in stable (case BB) mass transfer onto a stripped helium burning secondary.
3. The primary, now a naked metal star, explodes in a SN and becomes a NS.
4. There is a final highly non-conservative stable (case BB) mass transfer episode from the stripped post-helium-burning secondary onto the NS.
5. The secondary then explodes as an USSN, forming a DNS.
The two dominant channels, Channel I and Channel II, comprise 69% and 14% of all DNSs in our simulations (
$Z=0.0142$
), respectively. The remaining formation channels are mostly variations of the dominant channels. These variations either alter the sequence of events or avoid certain mass transfer phases. Some formation scenarios rely on fortuitous SN kicks. Some other exotic scenarios, which allow for the formation of DNS in which neither NS is recycled by accretion (e.g., Belczyński & Kalogera Reference Belczyński and Kalogera2001), comprise less than 2% of the DNS population.
3.2. Comparison with Vigna-Gómez et al. (Reference Vigna-Gómez2018)
This work generally uses similar assumptions and physics parameterisations as the preferred model of Vigna-Gómez et al. (Reference Vigna-Gómez2018), including the Fryer et al. (Reference Fryer, Belczynski, Wiktorowicz, Dominik, Kalogera and Holz2012) delayed SN engine. Although the qualitative results are similar, there are some quantitative changes due to updated model choices and corrections to the COMPAS population synthesis code as described in Section 2.2 and Appendix A. For example, the percentage of systems forming though Channel I remains 
${\approx}70 \%$
, but now only 
${\approx}14\%$
of systems experience Channel II, instead of 
${\approx}21 \%$
in Vigna-Gómez et al. (Reference Vigna-Gómez2018). The main change concerns the DNS rates which in this work are a factor of a few lower than those in the preferred model of Vigna-Gómez et al. (Reference Vigna-Gómez2018). Here, we estimate the formation rate of all (merging) DNS to be 85 (60) 
$\rm Gpc^{-3}/yr$
. Vigna-Gómez et al. (Reference Vigna-Gómez2018). reports the formation rate of all (merging) DNS to be 369(281) 
$\rm Gpc^{-3}/yr$
for the preferred model. We discuss rates in more detail in Section 4.4.6.
3.3. CEEs leading to DNS formation in the HR diagram
For all properties, we present a colour-coded HR diagram, normalised distribution, and cumulative distribution function (CDF). In Figure 2, we present our synthetic population of DNS progenitors at the onset of RLOF leading to a CEE. They are coloured according to the stellar type of the donor at RLOF, which is specified using the nomenclature from Hurley et al. (Reference Hurley, Pols and Tout2000)Footnote c. Additionally, Figure 2 shows the normalised distributions of luminosity (
$L_{\textrm{donor}}$
), effective temperature (
$T_{\textrm{eff,donor}}$
), and stellar type of the donor. In the case of a double-core CEE, the donor is defined as the more evolved star from the binary, which is the primary in Channel II.
Main properties of the donor star at the onset of RLOF leading to the CEE in DNS-forming binaries. Top: HR diagram coloured by stellar phase: HG (blue), GB (orange), CHeB (yellow), and EAGB (purple). The sizes of the markers represent their sampling weight. We show the progenitor of the luminous red nova M101 OT2015-1 (Blagorodnova et al. Reference Blagorodnova2017) with a star symbol. The solid black lines indicate ZAMS and TAMS loci for a grid of SSE models (Hurley et al. Reference Hurley, Pols and Tout2000) at 
$Z\approx0.0142$
. We show the evolution of a single non-rotating 
$16\ \textrm{M}_{\odot}$
star, from ZAMS to the end of the giant phase: the dotted dark grey line shows a MIST stellar track from Choi et al. (Reference Choi, Dotter, Conroy, Cantiello, Paxton and Johnson2016) and the dashed grey line shows the stellar track from Pols et al. (Reference Pols, Schröder, Hurley, Tout and Eggleton1998, Reference Pols, Schroeder, Hurley, Tout and Eggleton2009). The dash-dotted light blue and solid green lines show how fitting formulae from Hurley et al. (Reference Hurley, Pols and Tout2000) lead to a bifurcation after the MS for stars with masses between 12.9 and 13.0 
$M_\odot$
. This bifurcation is related to which stars are assumed to begin core helium burning while crossing of the HG or only after it: see the presence (lack) of the blue loop in the 
$12.9\ (13.0)\ \textrm{}M_{\odot}$
track. Grey lines indicate stellar radii of 
$R=\{10,100,500,1000\}\ \textrm{R}_{\odot}$
. Bottom: Normalised distributions in blue (left vertical axis) and CDF in orange (right vertical axis) of luminosity (left panel), effective temperature (middle panel) and stellar type (right panel). Black error bars indicate 
$1\sigma$
sampling uncertainty in the histograms. Grey lines show 100 bootstrapped distributions that indicate the sampling uncertainty in the CDFs. The CDFs show a subset of 365 randomly sampled values, which is the same number of DNS in our population, for each bootstrapped distribution.
In Figure 2, there is a visually striking feature: the almost complete absence of systems which forms a white polygon around 
$\log_{10}(L_{\textrm{donor}}/\textrm{L}_{\odot})=4.5$
. This feature is a consequence of the fitting formulae used for single stellar evolution (SSE) (c.f. Figures 14 and 15 of Hurley et al. Reference Hurley, Pols and Tout2000). This white region is bounded by the evolution of a 
$12.9\ \textrm{M}_{\odot}$
and a 
$13.0\ \textrm{M}_{\odot}$
star at 
$Z = 0.0142$
. The MS evolution of both stars is quite similar. After the end of the MS, there is a bifurcation point arising from the lower mass system experiencing a blue loop and the higher mass system avoiding it. This bifurcation is enhanced by the sharp change in the 
$T_{\textrm{eff}}-L$
slope from the interpolation adopted by Hurley et al. (Reference Hurley, Pols and Tout2000) during the HG phase. This change in slope around 
$\log_{10}(L_{\textrm{donor}}/\textrm{L}_{\odot})=4.5$
and 
$\log_{10}(T_{\textrm{eff,donor}}/\textrm{K})=4.4$
is model-dependent but we do expect to have some differences in the evolution of stars around that mass. This bifurcation corresponds to the transition around the First Giant Branch, which separates intermediate-mass and high-mass stars. Stellar tracks from Choi et al. (Reference Choi, Dotter, Conroy, Cantiello, Paxton and Johnson2016) also display a bifurcation point, but models and interpolation are smoother than those in Hurley et al. (Reference Hurley, Pols and Tout2000).
A rare example of how a system could end up in the forbidden region is the following. If a star experiences a blue loop, it contracts and then re-expands before continuing to evolve along the giant branch (GB). If the companion experiences a SN with a suitable kick while the star is in this phase and the orbit is modified appropriately in the process, the system may experience RLOF. A fortuitous kick making the orbit smaller, more eccentric, or both would be an unusual but not implausible outcome of a SN.
3.4. Properties of the donor
We report the luminosity, effective temperature, stellar phase, mass and core mass fraction of the donor (
$f_{\textrm{donor}}\equiv m_{\textrm{core,donor}}/m_{\textrm{donor}}$
), as presented in Table 1. The luminosity and effective temperature limits are 
$\log_{10}\ [L_{\textrm{donor,min}}/\textrm{L}_{\odot},L_{\textrm{donor,max}}/\textrm{L}_{\odot}] = [ \lowLimLuminosity,\uppLimLuminosity]$
and 
$\log_{10}\ [T_{\textrm{eff,donor,min}}/\textrm{K},T_{\textrm{eff,donor,max}}/\textrm{K}] = [\lowLimTemperature,\uppLimTemperature]$
, respectively. In Figure 2, we highlight the stellar phase, which is colour-coded. While the evolution in the HR diagram is itself an indicator of the evolutionary phase of the star, our stellar models follow closely the stellar-type nomenclature as defined in Hurley et al. (Reference Hurley, Pols and Tout2000). Donors which engage in a CEE leading to DNS formation can be in the HG (4%), GB (7%), CHeB (59%), or EAGB (30%) phase.
Properties of the donor star and the binary at the onset of RLOF leading to a CEE. In this Table, we list the symbols and units for each parameter, as well as the figure where the parameter is presented.
In the case of Channel I, donors are HG or CHeB stars; they span most of the parameter space from terminal-age MS (TAMS) until the end of CHeB, with a temperature range of 
$\sim 0.5$
dex. In the case of Channel II, donors are GB or EAGB giant-like stars. The parameter space in the HR diagram for these giant-like donors is significantly smaller, spanning an effective temperature range of only 
$\sim 0.1$
dex.
The limits in the mass of the donor are 
$[m_{\textrm{donor,min}},$
$ m_{\textrm{donor,max}}] =[\lowLimMass, \uppLimMass]\ \textrm{M}_{\odot}$
. The core mass fraction, shown in Figure 3, has limits of 
$[f_{\textrm{core,donor,min}},$
$ f_{\textrm{core,donor,max}}]=[\lowLimCoreMassFraction, \uppLimCoreMassFraction]$
. The core mass fraction can serve as a proxy for the evolutionary phase.
Pre-CEE donor properties of all DNS-forming systems: mass (top), core mass fraction (middle), and envelope binding energy (bottom). The core mass fraction is defined as 
$f_{\textrm{core,donor}}\equiv m_{\textrm{core,donor}}/m_{\textrm{donor}}$
. In the case of a double-core CEE, the binding energy is the sum of the individual envelope binding energies. Yellow systems with binding energies larger than 
$\log_{10}\ |E_{\textrm{bind}}/\textrm{erg}|\approx 48.5$
during the red supergiant phase are double-core CEE systems. For more details, see Section 3.4. See the caption of Figure 2 for further explanations.
Pre-CEE orbital properties of all DNS-forming systems. The binary properties presented are eccentricity (top) and semi-major axis (bottom). The orbital properties do not account for tidal circularisation. For more details, see Section 3.5. See the caption of Figure 2 for further explanations.
We report the binding energy of the envelope (see Figure 3) as defined in Equation 4. In the case of a double-core CEE, the binding energy of the CE is assumed to be 
$E_{\textrm{bind}}=E_{\textrm{bind,donor}}+E_{\textrm{bind,comp}}$
. The binding energy falls in the range 
$\log_{10}\ ({-}[E_{\textrm{bind,min}},E_{\textrm{bind,max}}]/\textrm{erg}) =[\uppLimEbind,\lowLimEbind]$
. The systems with the most tightly bound envelopes, and therefore the lowest (most negative) binding energies, are those experiencing a CEE shortly after TAMS or as double-core CEEs. For double-core systems, the envelope of the less evolved companion is more bound than the one of the more evolved donor star.
3.5. Properties of the binary
We also report the properties of each binary by colour-coding the property of interest in the HR diagram. We report the eccentricity, semi-major axis, total mass (
$m_{\textrm{total}} = m_{\textrm{comp}}+m_{\textrm{donor}}$
), mass ratio (
$q = m_{\textrm{comp}}/m_{\textrm{donor}}$
), and the ratio of the circularisation timescale (
$\tau_{\textrm{circ}}$
) to the radial expansion timescale (
$\tau_{\textrm{radial}}$
), as presented in Table 1. All quantities are reported at the onset of the RLOF unless stated otherwise.
The eccentricity, semi-major axis, and masses of the system determine the orbital energy and angular momentum of a binary (in the two point mass approximation). The eccentricity and semi-major axis distributions shown in Figure 4 do not account for tidal circularisation. The eccentricities span the entire allowed parameter range 
$0 \le e < 1$
. The eccentricity distribution has a sharp feature around 
$e\approx 0$
. Systems with 
$e\approx 0$
are typically those from Channel II, where the double-core CEE happens as the first mass transfer interaction, without any preceding SN to make the binary eccentric given our assumption of initially circular binaries (further discussion of this choice is in Section 4.4.3). Meanwhile, the most eccentric binaries have the smallest periapses and interact the earliest during the evolution of the donor, explaining the trend of greater eccentricities for smaller donor sizes in Figure 4.
The semi-major axis distribution, shown in Figure 4, has limits of 
$[a_{\textrm{min}},a_{\textrm{max}}]=[\lowLimSemiMajorAxis, \uppLimSemiMajorAxis]\ \textrm{R}_{\odot}$
. The very few extremely wide systems correspond to very eccentric binaries, almost unbound during the SN explosion (e.g., 
$e \approx 0.9999$
for the widest binary). While those limits are broad, the limits in periastron are 
$[a_{\textrm{p,min}},a_{\textrm{p,max}}] \approx [\lowLimPeriastron,\uppLimPeriastron]\ \textrm{R}_{\odot}$
. (Very rarely, even smaller periapses are possible when fortuitous SN kicks send the newly formed NS plunging into the envelope of an evolved companion on a very eccentric orbit; however, it is not clear whether such events lead to a CEE or to a more exotic outcome, such as the formation of a Thorne & Żytkow Reference Thorne and Żytkow1977 object). The total mass distribution, shown in Figure 5, has limits of 
$[m_{\textrm{total,min}},m_{\textrm{total,max}}]=[\lowLimTotalMass, \uppLimTotalMass]\ \textrm{M}_{\odot}$
.
Pre-CEE mass of all DNS-forming systems. The binary properties presented are total mass (top) and mass ratio (bottom). For more details, see Section 3.5. See the caption of Figure 2 for further explanations.
We compute the mass ratio at the onset of the RLOF leading to the CEE. The mass ratio, shown in Figure 5, has limits of 
$[q_{\textrm{min}},q_{\textrm{max}}]=[\lowLimMassRatio, \uppLimMassRatio]$
. The broad distribution in fact consists of two distinct peaks, one close to 
$q=0$
and the other close to 
$q=1$
, with a large gap between 
$0.18 \leq q \leq 0.98$
(see Figure 5). The extreme mass ratio systems correspond to CEEs from Channel I, where the companion is a NS. The 
$q \approx 1$
systems correspond to CEEs from Channel II, where there is a double-core CEE with a non-compact companion star. The systems with 
$q > 1$
are double-core CE systems with 
$q_{\textrm{ZAMS}} \approx 1$
which, at high metallicity, may reverse their mass ratio via mass loss through winds before the primary star expands and undergoes RLOF.
3.6. Tidal timescales in pre-CE systems
Given the uncertainties in the treatment of tides, and our interest in comparing the impact of different tidal prescriptions as discussed below (see Section 4.4.2), we do not include tidal synchronisation or circularisation in binary evolution modelling for this study. Instead, we consider whether tides would be able to efficiently circularise the binary before the onset of RLOF leading to a CEE. As discussed in Section 2.3.3, we use the ratio of the circularisation timescale to the radial expansion timescale as a proxy for the efficiency of tidal circularisation of an expanding star about to come into contact with its companion. If 
$\tau_{\textrm{circ}}/\tau_{\textrm{radial,donor}} > 1$
, we label the binary as still eccentric at RLOF. Given that the circularisation timescale is longer than the synchronisation timescale (see Section 2.3), we focus on the former and assume that if the binary is able to circularise, it will already be synchronous.
Figure 6 shows the ratio 
$\tau_{\textrm{circ}}/\tau_{\textrm{radial,donor}}$
under our default assumption in which both HG and CHeB stars have fully convective envelopes for the purpose of tidal circularisation calculations and experience the equilibrium tide. This assumption results in 82% of the systems being circular at the onset of RLOF.
The analysis of circularisation timescales is mostly relevant for systems formed through Channel I (see Section 4.4.2). There are two reasons for this. The first one is that they are expected to acquire a non-zero eccentricity after the first SN. The second one is that they are more likely to have a radiative or only partially convective envelope, making circularisation less efficient. On the other hand, systems formed through Channel II have a fully convective envelope, which allows for efficient tidal circularisation and synchronisation. We apply the low-eccentricity approximation described in Section 2.3.1 to computing the tidal timescales of these systems, even though they are circular by construction before the first SN (we discuss this assumption in Section 4.4.3).
4. Discussion
In this section, we discuss the properties of CEEs experienced by isolated stellar binaries evolving into DNSs and present some of the caveats in our COMPAS rapid population synthesis models.
4.1. CEE subpopulations in evolving DNSs
4.1.1. Formation channels
There are two main formation channels leading to DNS formation. Channel I involves high-mass ratio single-core CEE between a NS primary and a post-MS secondary. Channel I has been studied thoroughly in the literature, for example, Bhattacharya & van den Heuvel (Reference Bhattacharya and van den Heuvel1991); Tauris & van den Heuvel (Reference Tauris and van den Heuvel2006); Tauris et al. (Reference Tauris2017) and references therein. Channel II involves a double-core CE between two post-MS stars. A similar channel has been proposed by Brown (Reference Brown1995) and Dewi et al. (Reference Dewi, Podsiadlowski and Sena2006), among others. Channel II requires similar masses at ZAMS driven by the need of similar evolutionary timescales so that both stars are post-MS giants at the time of their first interaction. For low- (high-) mass stars, the difference in ZAMS mass can be up to 3 (7)%, in agreement with Dewi et al. (Reference Dewi, Podsiadlowski and Sena2006). Our Channel II has an additional case BB mass transfer episode from a helium shell burning primary onto a helium MS secondary.
4.1.2. Subpopulations and tidal circularisation
We separate CEE donors into three distinct subpopulations depending on their evolutionary phase at the onset of RLOF: giants, cool, and hot (see Table 2 and Figure 7).
Distinct DNS subpopulations as described in Section 4.1 and presented in Figure 7.
Ratio of tidal circularisation timescale to the star’s radial expansion timescale for all DNS-forming systems. We present the default scenario where all evolved stars, including HG and CHeB stars, are assumed to have formed a fully convective envelope. If 
$\log_{10}(\tau_{\textrm{circ}}/\tau_{\textrm{radial}})\le 0$
, we assume that binaries circularise before the onset of the CEE. Binaries indicated with blue (red) dots are predicted to have circular (eccentric) orbits. We cap 
$\ -2 \le \log_{10}(\tau_{\textrm{circ}}/\tau_{\textrm{radial}}) \le 2$
to improve the plot appearance. The grey shaded region in the histogram highlights the systems which circularise by the onset of RLOF. For more details, see Section 3.6. See the caption of Figure 2 for further explanations.
DNS-forming binaries clustered by the donor type at the onset of the CEE. Subpopulations: (a) giant donors with fully convective envelopes in blue, (b) HG or CHeB donors with partially convective envelopes in red, and (c) HG or CHeB donors which have not yet formed a deep convective envelope in yellow. For more details, see Section 3.4. See the caption of Figure 2 for further explanations.
The first one, giants, correspond to giant donors with fully convective envelopes. The other two subpopulations correspond to HG or CHeB donors, most of them evolving via the single-core Channel I. We distinguish between cool donors with a partially convective envelope and hot donors with a radiative envelope. We follow Belczynski et al. (Reference Belczynski2008) in using the temperature 
$\log_{10}\ (T_{\textrm{eff,donor}}/\textrm{K}) = 3.73$
as the boundary between the cool and hot subpopulations.
The presence and depth of a convective envelope impacts the response of the star to mass loss and, hence, the dynamical stability of mass transfer. In particular, hot donors lacking a deep convective envelope may be stable to mass transfer and avoid a CEE. At the same time, some of the less evolved hot donors may not survive a CEE even if they do experience dynamical instability (pessimistic variation). Klencki et al. (Reference Klencki, Nelemans, Istrate and Chruslinska2020) use detailed stellar evolution models to argue that only red supergiant donors with deep convective envelopes are able to engage in and survive a CEE. For their assumptions, this would reduce the estimated rate of DNS formation. However, and similar to this study, they focus on RLOF structures and not the structures at the moment of the instability.
Here, we focus only on the impact of the assumed structure of the donor on the efficiency of tidal circularisation and do not account for possible consequences for mass transfer stability. We compare three alternative models in Figure 8.
CDF of the ratio of the circularisation timescale to the donor radial expansion timescale computed at RLOF onset leading to CEE for all DNS-forming systems. Here, we present three scenarios. The solid blue line is our default assumption: all donors have a deep convective envelope (same as in left panel of Figure 6). The red dashed line follows Hurley et al. (Reference Hurley, Tout and Pols2002) with the assumption that CHeB tidal evolution is dominated by the dynamical tide, that is, that CHeB stars have a radiative envelope. The yellow dotted line follows Belczynski et al. (Reference Belczynski2008) in assuming that stars with 
$\log\ T_{\textrm{eff}}\le 3.73\ \textrm{K}$
have a fully convective envelope, for both HG and CHeB donors; and a fully radiative envelope otherwise, as in Figure 7. For more details, see Section 3.6.
All DNS-forming binaries from our Fiducial model are shown here. We present the post-CEE separation 
$a_{\textrm{f}}$
as a function of the absolute value of the envelope binding energy 
$|E_{\textrm{bind}}|$
. For the double-core scenario, the binding energy is 
$E_{\textrm{bind}}=E_{\textrm{bind,donor}}+E_{\textrm{bind,comp}}$
. The size of the marker indicates the sampling weight and its colour shows the mass ratio q. This Figure can be compared to Figures 1 and 2 from Iaconi & De Marco (2019). That study presents simulations of CE binaries and observations of post-CE binaries. Most systems presented here do not feature in Iaconi & De Marco (2019).
Our default tidal circularisation model assumes that all evolved donors, including both HG and CHeB stars, have fully convective envelopes and therefore experience efficient equilibrium tides. Our default assumption estimates that 
${\rm{18}}$
% of systems will be eccentric at the onset of the RLOF leading to the CEE. This is the lowest fraction of eccentric systems among all variations because tides are particularly efficient for stars with convective envelopes.
In reality, CHeB stars are expected to begin the CHeB phase with a radiative envelope and develop a deep convective envelope by the end of it. The single stellar fits from Hurley et al. (Reference Hurley, Pols and Tout2000) do not contain explicit information about the moment when this transition occurs. Hurley et al. (Reference Hurley, Tout and Pols2002) assume that all CHeB stars have a radiative envelope and that the dynamical tide is dominant in their tidal evolution. Adopting this assumption leads to 
${\rm{68}}$
% of binaries remaining eccentric at the onset of the RLOF leading to the CEE.
Alternatively, Belczynski et al. (Reference Belczynski2008) assume that hot stars with 
$\log_{10}\ (T_{\textrm{eff}}/\textrm{K}) > 3.73$
have a radiative envelope, while cool stars with 
$\log_{10}\ (T_{\textrm{eff}}/\textrm{K}) \le 3.73$
have a convective envelope. Adopting this assumption leads to 
${\rm{40}}$
% of binaries remaining eccentric at the onset of the RLOF leading to the CEE.
According to our estimates, a significant fraction of systems will be eccentric at RLOF. These estimates were made within the framework of the fitting formulae for SSE from Hurley et al. (Reference Hurley, Pols and Tout2000). More detailed fitting formulae, which include the evolutionary stage of stars as well as the mass and radial coordinates of their convective envelopes, would allow for a self-consistent determination of whether a star has a radiative, a partially convective, or a fully convective envelope for both dynamical stability and tidal circularisation calculations.
4.2. CEEs as candidates for luminous red novae transients
Recently, the luminous red nova transient M101 OT2015-1 was reported by Blagorodnova et al. (Reference Blagorodnova2017). This event is similar to other luminous red novae associated with CEEs (Ivanova et al. Reference Ivanova, Justham, Avendano Nandez and Lombardi2013b). Following the discovery of M101 OT2015-1, archival photometric data from earlier epochs were found. Blagorodnova et al. (Reference Blagorodnova2017) used these to derive the characteristics of the progenitor. The inferred properties of the progenitor of M101 OT2015-1 are a luminosity of 
$L_{\textrm{donor}}\approx 87\,000\ \textrm{L}_{\odot}$
, an effective temperature of 
$T_{\textrm{eff,donor}}\approx 7\,000\ \textrm{K}$
and a mass of 
$m_{\textrm{donor}}=18 \pm 1\ \textrm{M}_{\odot}$
(see Figure 2 for location in the HR diagram).
Binary separations at CEEs leading to DNSs at the onset of RLOF (left) and after the CEE (right). We show the donor (
$m_{\textrm{donor}}$
) and companion (
$m_{\textrm{comp}}$
) mass in both plots, with a solid grey line indicating 
$m_{\textrm{donor}}=m_{\textrm{comp}}$
. The colour bars, with different scales, show the pre-CEE periastron (left) and final semi-major axis (right).
Blagorodnova et al. (Reference Blagorodnova2017) found that the immediate pre-outburst progenitor of M101 OT2015-1 was consistent with an F-type yellow supergiant crossing the HG. If we take the inferred values for this star as the values at the onset of RLOF, then this star is consistent with pre-CEE stars in our predicted distribution of DNS-forming systems. However, we emphasise that the appearance of the donor star can change significantly between the onset of RLOF, that is, the point at which the models shown in Figure 2 are plotted, and dynamical instability.
Howitt et al. (Reference Howitt2020) explored population synthesis models of luminous red novae. Here, we use the same pipeline adopted for that study to explore the connection to DNS populations. Doing so, fewer than 0.02% of all luminous red novae lead to DNSs. These are among the most energetic luminous red novae and would be overrepresented in the magnitude-limited observable population. Future DNSs constitute nearly 10% of the subpopulation of luminous red novae with predicted plateau luminosities greater than 
$10^7\ \textrm{L}_{\odot}$
.
4.3. Eccentric RLOF leading to a CEE
We predict that the subpopulation of giant donors with fully convective envelopes and cool donors with partially convective envelopes are likely to be circular at the onset of the CEE (see Figures 6 and 8). On the other hand, we find that the subpopulation of hot donors often does not circularise by the onset of the CEE. This subpopulation with hot donors are binaries with high eccentricities at the onset of the RLOF (see Figure 4).
This result raises questions about the initial conditions of a CEE, which is often assumed to begin in a circular orbit, both in population synthesis studies and in detailed simulations. Population synthesis codes such as SEBA (Portegies Zwart & Verbunt Reference Portegies Zwart and Verbunt1996; Portegies Zwart & Yungelson Reference Portegies Zwart and Yungelson1998; Toonen, Nelemans, & Portegies Zwart Reference Toonen, Nelemans and Portegies Zwart2012), STARTRACK (Belczynski et al. Reference Belczynski, Kalogera and Bulik2002, Reference Belczynski2008), binary stellar evolution (BSE) (Hurley et al. Reference Hurley, Tout and Pols2002), the Brussels code (De Donder & Vanbeveren Reference De Donder and Vanbeveren2004), COMPAS (Stevenson et al. Reference Stevenson, Vigna-Gómez, Mandel, Barrett, Neijssel, Perkins and de Mink2017; Vigna-Gómez et al. Reference Vigna-Gómez2018), ComBinE (Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018), and customised software based on them all assume that RLOF commences in circular binaries. Detailed simulations, such as those of Passy et al. (Reference Passy2012), MacLeod, Ostriker, & Stone (Reference MacLeod, Ostriker and Stone2018), and others, often make the assumption of an initially circular orbit (but see Staff et al. Reference Staff, De Marco, Macdonald, Galaviz, Passy, Iaconi and Mac Low2016, discussed in more detail in Section 4.3.2).
4.3.1. Theory of mass transfer in eccentric binaries
Mass transfer in eccentric binaries has been explored with both semi-analytical and analytical methods (Matese & Whitmire Reference Matese and Whitmire1983, Reference Matese and Whitmire1984; Sepinsky et al. Reference Sepinsky, Willems, Kalogera and Rasio2007, Reference Sepinsky, Willems, Kalogera and Rasio2009, Reference Sepinsky, Willems, Kalogera and Rasio2010; Dosopoulou & Kalogera Reference Dosopoulou and Kalogera2016a,b).
The analysis of Sepinsky et al. (Reference Sepinsky, Willems, Kalogera and Rasio2007) et al. assumes fully conservative mass transfer. While they consider mass transfer from a stellar donor onto a NS, this assumption and the typical 
$10^{-9}\ \textrm{M}_{\odot}/yr$
mass transfer rate they consider is relevant for low-mass X-ray binaries, not DNS progenitors as discussed here.
Dosopoulou & Kalogera (Reference Dosopoulou and Kalogera2016b) study orbital evolution considering both conservative and non-conservative mass transfer. The latter scenario is particularly relevant for DNS formation. We assume a mass transfer rate of 
$10^{-5}\ \rm M_{\odot}/yr$
, a 
$1.44\ \rm M_{\odot}$
NS and mass loss from the vicinity of the NS (isotropic re-emission) as parameters in their Equation (44). Highly eccentric systems (
$e>0.9$
) have circularisation timescales of more than 1 Myr due to mass transfer, under their assumption that all mass transfer happens at periapsis. This timescale is reduced to around a 1000 yr for mass transfer in 
$e<0.1$
binaries, although this can be very sensitive to assumptions about the specific angular momentum lost at the level of an order of magnitude. The assumption of instantaneous mass transfer at periapsis is questionable for low-eccentricity binaries, precisely those which may efficiently circularise through mass transfer.
Hamers & Dosopoulou (Reference Hamers and Dosopoulou2019) noted that evolution towards circularisation from Dosopoulou & Kalogera (Reference Dosopoulou and Kalogera2016b) could lead to (non-physical) negative eccentricity solutions. They proposed a revised analytic model for mass transfer in eccentric binaries. This study takes into account the separation and eccentricity evolution of an initially eccentric system at RLOF. However, their model is only valid in the regime of fully conservative mass transfer and is therefore more restricted than the general formalism Dosopoulou & Kalogera (Reference Dosopoulou and Kalogera2016b). Mass transfer episodes in binaries which will become DNSs are typically non-conservative. Mass transfer from a post-MS donor onto a MS companion, such as the first mass transfer episode from Channel I, is generally only partly conservative (Schneider et al. Reference Schneider, Izzard, Langer and de Mink2015). Mass transfer onto a NS companion is highly non-conservative, almost in the fully non-conservative limit (Tauris et al. Reference Tauris, Langer and Podsiadlowski2015).
A full understanding of the evolution of eccentric systems in RLOF is yet to be achieved. A detailed treatment of non-conservative mass transfer in an eccentric binary could yield different criteria for dynamical stability and, ultimately, for determining if a system engages in a CEE.
4.3.2. Modelling of mass transfer in eccentric binaries
Numerical methods and simulations have also been used to study mass transfer in eccentric binaries (Bobrick, Davies, & Church Reference Bobrick, Davies and Church2017; Church et al. Reference Church, Dischler, Davies, Tout, Adams and Beer2009; Lajoie & Sills Reference Lajoie and Sills2011; Regös, Bailey, & Mardling Reference Regös, Bailey and Mardling2005; Staff et al. Reference Staff, De Marco, Macdonald, Galaviz, Passy, Iaconi and Mac Low2016; van der Helm, Portegies Zwart, & Pols Reference van der Helm, Portegies Zwart and Pols2016).
Staff et al. (Reference Staff, De Marco, Macdonald, Galaviz, Passy, Iaconi and Mac Low2016) carried out hydrodynamic simulations of a 
${\approx} 3\ \textrm{M}_{\odot}$
giant star with a less massive MS companion in an eccentric orbit. They conclude that eccentric systems transfer mass only during the periastron passage, which delays the onset of the CEE. Each periastron passage also makes the binary less eccentric.
Gilkis, Soker, & Kashi (Reference Gilkis, Soker and Kashi2019) discuss the passage of a NS through the envelope of a giant star, likely on an eccentric orbit, and conclude that the system might be able to eject the envelope or lead to a merger between the NS and the core, which they call a CE jets SN (Soker & Gilkis Reference Soker and Gilkis2018, see also Schrøder et al. 2020). The former results in a less luminous transient, comparable in energetic to luminous red novae (Kashi & Soker Reference Kashi and Soker2016). The interaction is driven by jets which might enhance mass loss at periastron passages, keeping the system eccentric (Kashi & Soker Reference Kashi and Soker2018). These jets might, in some cases, prevent CEEs (Shiber & Soker Reference Shiber and Soker2018).
4.3.3. Observations of eccentric mass-transferring binaries
Eccentric semi-detached and contact binaries, that is, mass-transferring binaries, have been previously observed in low-mass systems (Petrova & Orlov Reference Petrova and Orlov1999). Eccentric (
$e\lesssim 0.2$
) MS White Dwarf binaries which are believed to have experienced RLOF are not rare (Kawahara et al. Reference Kawahara, Masuda, MacLeod, Latham, Bieryla and Benomar2018; Masuda et al. Reference Masuda, Kawahara, Latham, Bieryla, Kunitomo, MacLeod and Aoki2019; Vos et al. Reference Vos, Østensen, Németh, Green, Heber and Van Winckel2013). Jayasinghe et al. (Reference Jayasinghe, Stanek, Kochanek, Thompson, Shappee and Fausnaugh2019) found a more massive B-type Heartbeat star in an eccentric (
$e=0.58$
) orbit. Heartbeat stars exhibit clear signatures of tidal oscillations at each periastron passage. While there is no evidence for accretion, it is likely that the system reported in Jayasinghe et al. (Reference Jayasinghe, Stanek, Kochanek, Thompson, Shappee and Fausnaugh2019) will engage in RLOF at some later point.
Sirius (Gatewood & Gatewood Reference Gatewood and Gatewood1978; van den Bos Reference van den Bos1960) is a MS White Dwarf binary with 
$e=0.59$
which, according to canonical binary evolution dynamics, should have circularised when the White Dwarf progenitor became a giant. Bonacić Marinović, Glebbeek, & Pols (Reference Bonačić Marinović, Glebbeek and Pols2008) propose a model which allows for tides, mass loss, and mass transfer in an eccentric orbit, physically motivated by Sirius. Following Bonacić Marinović et al. (Reference Bonačić Marinović, Glebbeek and Pols2008), Saladino & Pols (Reference Saladino and Pols2019) carried out hydrodynamic simulations of binary stars with significant wind-driven mass loss and find that this eccentricity-enhancing mechanism is non-negligible.
4.3.4.
$\gamma ^2$
Velorum as an eccentric massive post-mass-transferring binary
North et al. (Reference North, Tuthill, Tango and Davis2007) reported the orbital solution and fundamental parameter determination of the massive binary 
$\gamma^2$
Velorum. This binary has a reported period of 
$78.53\,\pm\,0.01$
d and an eccentricity of 
$0.334\,\pm\,0.003$
, with an inferred mass of 
$28.5\,\pm\,1.1\ \textrm{M}_{\odot}$
for the O-star primary and 
$9.0\pm0.6\ \textrm{M}_{\odot}$
for the Wolf–Rayet secondary. While 
$\gamma^2$
Velorum did not experience a CEE, it could have experienced some mass transfer as an eccentric system. Eldridge (Reference Eldridge2009) discusses 
$\gamma^2$
Velorum as post-mass transfer binary system. In that work, Eldridge (Reference Eldridge2009) takes into account how the evolutionary stage of the donor during mass transfer determines the efficiency of tidal circularisation. They point out that a less evolved star with a radiative envelope is not likely to circularise during the mass transfer phase. This would lead to a post-mass transfer eccentric system such as 
$\gamma^2$
Velorum.
4.4. Caveats and limitations
4.4.1. CEE and delayed dynamical instability
The uncertainties in our stellar and binary models propagate to uncertainties in whether a mass-transferring system experiences a CEE. We compare the response of the radius of the donor to mass loss to the response of the orbit to mass transfer to determine whether a binary experiences a CEE (see Section 2.2). This approach relies on determining the appropriate response of the donor to (adiabatic) mass loss, the amount of mass that the companion can accrete, and the specific angular momentum removed from the binary by the non-accreted mass; all of these quantities have uncertainties and are model-dependent. Other population synthesis codes directly use the mass ratio at RLOF to determine whether the mass transfer will be stable (e.g., Claeys et al. Reference Claeys, Pols, Izzard, Vink and Verbunt2014; Hurley et al. Reference Hurley, Tout and Pols2002).
The evolution of a mass-transferring system is non-trivial. One possibility is delayed dynamical instability, in which the donor experiences a prolonged mass transfer phase before it becomes dynamically unstable (Ge et al. Reference Ge, Hjellming, Webbink, Chen and Han2010; Hjellming & Webbink Reference Hjellming and Webbink1987; Ivanova & Taam Reference Ivanova and Taam2004). This can lead the donor to be significantly under-luminous at the moment when the CEE begins, compared to its appearance at the onset of the mass transfer episode itself (Podsiadlowski, Rappaport, & Pfahl Reference Podsiadlowski, Rappaport and Pfahl2002). We report the properties at the onset of the RLOF because we do not account for delayed dynamical instability.
The opposite is also possible, where initially unstable systems may reach a stable configuration after ejecting only a fraction of the CE. Pavlovskii et al. (Reference Pavlovskii, Ivanova, Belczynski and Van2017) found that some massive giant donors with stellar mass BH companions, which were previously expected to experience a CEE, might experience stable mass transfer instead.
In general, the transition from stable to unstable mass transfer is not fully understood. The details of the threshold of stability are regulated by the hydrodynamics of the material lost from the donor and the thermal response of the donor. This confluence of physical processes will always be at play because, at the low-mass exchange rates that mass loss is initiated at, both dynamical and thermal readjustments of the binary are taking place on competing timescales. Accessing this regime in three-dimensional hydrodynamic simulations is challenging, as discussed by Reichardt et al. (Reference Reichardt, De Marco, Iaconi, Tout and Price2019) because near-stable configurations evolve over long timescales and are susceptible to small numerical perturbations. In the case of eccentric orbits, dynamical tides raised on the donor also play a significant role and the eventual orbital evolution is determined by a combination of tidal, thermal, and hydrodynamic evolution.
4.4.2. Tidal evolution
Both the tidal circularisation timescales and the applicability of various types of tides are highly uncertain. As discussed in Sections 2.3 and 4.1, we explore the impact of assumptions about the dominant tidal mechanism based on the evolutionary phase of the donor by considering the tidal circularisation timescale at the onset of RLOF. While this approach makes it possible to analyse the impact of different choices without re-analysing the full population, it does mean that tides are not self-consistently included throughout the evolution of the binary.
We also make a number of simplifying assumptions about the efficiency of tidal circularisation. For example, we apply the equilibrium tide to convective envelope donors regardless of the orbital eccentricity, although the perturbation-from-equilibrium approximation is unlikely to be valid for very eccentric binaries which are not pseudo-synchronised at periapsis. We crudely approximate coefficients in the tidal circularisation timescale equations based on Rasio et al. (Reference Rasio, Tout, Lubow and Livio1996), Hurley et al. (Reference Hurley, Tout and Pols2002), and references therein.
Witte & Savonije (Reference Witte and Savonije1999a,b) discuss how resonance locking could enhance pre-RLOF circularisation of a 
$10\ \textrm{M}_{\odot}$
MS star with a 
$1.4\ \textrm{M}_{\odot}$
NS companion. This system is similar to the phase immediately after the first SN in Channel I. The timescales on which resonance locking occurs are typically a few million years and could lead to less eccentric orbits at the onset of RLOF. It is uncertain how much of an effect resonance locking would have on a population of massive interacting binaries.
4.4.3. Zero-eccentricity initial distribution
We assume that binaries are circular at birth. This assumption is justified for close binaries, which are tidally circularised at birth, but is not consistent with observations of wide binaries (see Abt, Gomez, & Levy Reference Abt, Gomez and Levy1990; Levato et al. Reference Levato, Malaroda, Morrell and Solivella1987; Kobulnicky et al. Reference Kobulnicky2014; Sana et al. Reference Sana2012 as presented in Figure 3 of Moe & Di Stefano Reference Moe and Di Stefano2017). Our goal is to be conservative when studying eccentricity at the onset of the CEE. Thus all changes in eccentricity from an initially circular binary are due to the subsequent binary evolution. Vigna-Gómez et al. (Reference Vigna-Gómez2018) showed that using a thermal eccentricity distribution at birth decreases the DNS formation rate by about a half but has no significant effect on the orbital properties of DNS (we follow the eccentricity distribution from Vigna-Gómez et al. Reference Vigna-Gómez2018 to be able to make a more direct comparison with those results). However, the eccentricity distribution at ZAMS likely affects the eccentricity distribution at the onset of RLOF. This is particularly true for Channel II binaries, which enter the double-core CEE without previous interactions or SNe, thus retaining their birth eccentricity modulo tidal effects. The current 
$e=0$
peak associated with this system (see Figure 4) would be replaced by the birth eccentricity distribution.
4.4.4. Massive binary stars
In this work, we focused on CEEs during the formation of DNSs. Similar evolutionary pathways are experienced by other massive stellar binaries, including progenitors of BH–NS or BH–BH systems (Dominik et al. Reference Dominik2012; Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018; Neijssel et al. Reference Neijssel2019). The impact of eccentric RLOF is not in the exclusive interest of DCO formation. The role of tidal evolution and dynamical instability is fundamental for massive stellar mergers (Justham, Podsiadlowski, & Vink Reference Justham, Podsiadlowski and Vink2014; Podsiadlowski, Joss, & Hsu Reference Podsiadlowski, Joss and Hsu1992; Vigna-Gómez et al. Reference Vigna-Gómez, Justham, Mandel, de Mink and Podsiadlowski2019).
4.4.5. Dynamics
We do not consider the impact of dynamical interactions on the formation of DNSs. This could take the form of dynamically induced mergers in dense stellar environments, such as globular clusters (Andrews & Mandel Reference Andrews and Mandel2019 but see, e.g., Ye et al. Reference Ye, Fong, Kremer, Rodriguez, Chatterjee, Fragione and Rasio2020). Meanwhile, Kozai–Lidov oscillations (Kozai Reference Kozai1962; Lidov Reference Lidov1962) in hierarchical triple systems can drive up the eccentricity of the inner binary (see Naoz et al. Reference Naoz, Fragos, Geller, Stephan and Rasio2016 for a review) and contribute to the formation of merging DNSs (Hamers & Thompson Reference Hamers and Thompson2019). Both types of dynamical encounters can change the binary orbital evolution, including the eccentricity.
4.4.6. DNS merger rates
The merger rate of DNSs was inferred to fall in the range 110–3840 Gpc
$^3$
/yr with 90% confidence (Abbott et al. Reference Abbott2019) based on a single detection of GW170817 with a flat-in-rate prior (Abbott et al. Reference Abbott2017). The detection of GW190425 under the assumption of a DNS progenitor updates the local DNS merger rate to 250–2810 Gpc
$^3$
/yr (Abbott et al. Reference Abbott2020).
Vigna-Gómez et al. (Reference Vigna-Gómez2018) predicted a DNS merger rate of 
${\approx}280$
Gpc−3/yr. After several changes in COMPAS (see Section 3.2), the DNS merger rate is now 
${\approx}60$
Gpc−3/yr. This is more in line with the reported rates from other rapid population synthesis codes (Giacobbo, Mapelli, & Spera Reference Giacobbo, Mapelli and Spera2018; Giacobbo & Mapelli Reference Giacobbo and Mapelli2019b; Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018). Chruslinska et al. (Reference Chruslinska, Belczynski, Klencki and Benacquista2018) pointed out that all population synthesis models struggle to jointly predict the binary BH and DNS merger rate, particularly under standard assumptions. Most rapid population synthesis studies agree that predicted rates are difficult to reconcile with the DNS merger rates inferred from GW observations, particularly if future observations push these above 
$\gtrsim 1000\ \rm Gpc^{-3}/yr$
. Some particular choices of physics increase rates up to a few hundreds of DNSs 
$\rm Gpc^{-3}/yr$
, but the typical prediction is around a few tens of DNSs 
$\rm Gpc^{-3}/yr$
. Recently, Giacobbo & Mapelli (Reference Giacobbo and Mapelli2020) proposed that a revised natal kick prescription might prove important to reconcile DCO merger rates, as suggested by Chruslinska et al. (Reference Chruslinska, Belczynski, Klencki and Benacquista2018) (see more about natal kick prescriptions in Section 4.5).
The impact of the power law exponent (
$\alpha_{\rm IMF}$
) of the initial mass function on the DNS merger rate has been previously discussed by de Mink & Belczynski (Reference de Mink and Belczynski2015). They find that the rate can increase (decrease) by a factor of a few for plausible shallower (steeper) initial mass functions. Schneider et al. (Reference Schneider2018a, see also Farr & Mandel Reference Farr and Mandel2018; Schneider et al. Reference Schneider2018b) find that the initial mass function of young massive stars in the 30 Doradus region of the Large Magellanic Cloud may be shallower than the canonical Salpeter (Reference Salpeter1955) value. While this level of fluctuation in the initial mass function may prove insufficient to resolve the DNS rate discrepancy (Belczynski et al. Reference Belczynski2018), it may be one of the ingredients.
Observations of Galactic DNSs as radio pulsars point to a DNS merger rate that peaks at a few tens of mergers per Myr in the Galaxy, for example, 
$37_{-11}^{+24}$
per Myr according to Pol, McLaughlin, & Lorimer (Reference Pol, McLaughlin and Lorimer2020). This can be converted to a volumetric rate in units of Gpc−3/yr by multiplying the rate per Myr per Milky Way equivalent Galaxy by a factor of 
$\sim 10$
(Abadie et al. Reference Abadie2010). The peak of the rate extrapolated from Galactic observations thus falls between the typical binary population synthesis predictions and the peak of the rate inferred from GW observations. However, this extrapolation does not account for differences between the Galaxy and other environments (Abadie et al. Reference Abadie2010), and all rate intervals are broad.
Lau et al. (Reference Lau, Mandel, Vigna-Gómez, Neijssel, Stevenson and Sesana2020) use a synthesised DNS population to predict that a 4-yr Laser Interferometer Space Antenna (LISA) mission (Amaro-Seoane et al. Reference Amaro-Seoane2017; Baker et al. Reference Baker2019) will detect 35 Galactic DNSs. Andrews et al. (Reference Andrews, Breivik, Pankow, D’Orazio and Safarzadeh2020) predict between 46 and 240 Galactic DNSs for the same mission, depending on the assumed physical assumptions. LISA DNS observations will further constrain the DNS formation and merger rates.
4.4.7. Be X-ray binaries
The Be X-ray binary phase occurs in binaries consisting of a rapidly rotating Be star, a MS B star with emission lines from a decretion disc, and a NS which accretes from this disc. This phase can be one of the intermediate stages between ZAMS and DNS formation. We expect that most systems formed through Channel I experienced a Be X-ray binary configuration after the primary initiated the first mass transfer episode, spinning up the secondary, and exploded as a SN (see Figure 1).
The Small Magellanic Cloud (SMC) Be X-ray binary catalogue by Coe & Kirk (Reference Coe and Kirk2015) presents 69 systems. From those 69, 44 have an observed orbital period (
$P_{\textrm{orb}}<520$
d). Vinciguerra et al. (Reference Vinciguerra2020) used COMPAS to study Be X-ray binaries at SMC metallicity (
$Z\approx 0.0035$
). We focus on two particular variations from that study. Their default model follows similar COMPAS settings as the ones in this paper, while the preferred model allows for mass transfer to be more conservative than the default model, that is, half of the mass transferred is accreted by the companion during the mass transfer episode. The default model predicts 
$190\pm20$
Be X-ray binaries in the SMC; in this model, 
$\sim 59\%\ (46\%)$
of all (merging) DNS experience a Be X-ray binary phase. The preferred model predicts 
$80\pm10$
Be X-ray binaries in the SMC, with 
$\sim 96\%\ (98\%)$
of all (merging) DNS experiencing a Be X-ray binary phase. Both variations have a formation rate of merging DNSs of 
$\sim 10^{-5}\ \textrm{per\,M}_{\odot}$
. Vinciguerra et al. (Reference Vinciguerra2020) show that the preferred model better represents observations of Be X-ray binaries in comparison to the default model. However, there is no significant difference in the predicted DNS merger rate or the DNS mass distribution between these two models.
4.5. Comparison with other rapid population synthesis studies
It is challenging to compare results between population synthesis studies. The results depend on initial conditions, physical parameterisations, and computational methods. Most population synthesis software differ in the implementation of at least one of these, and often in all. Some of the arguably most important physical interactions, including SNe and mass transfer events, are treated differently between various research groups, codes, and particular projects.
In the context of SNe, there are several common choices for remnant mass and natal kick prescriptions. StarTrack, MOBSE, and COMPAS usually follow Fryer et al. (Reference Fryer, Belczynski, Wiktorowicz, Dominik, Kalogera and Holz2012) for the remnant mass prescription of compact objects, using either the rapid or delayed variation, while ComBinE use their own (Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018). The natal kick distribution varies not only between codes but also between papers. COMPAS initially followed StarTrack (e.g., Belczynski et al. Reference Belczynski2008, Reference Belczynski2018) with a no-kick model for electron-capture SN (ECSN) (Stevenson et al. Reference Stevenson, Vigna-Gómez, Mandel, Barrett, Neijssel, Perkins and de Mink2017) and more recently changed to a low-kick distribution (e.g., Vigna-Góómez et al. 2018). MOBSE has explored the impact of variations in the core-collapse SN (CCSN) and ECSN natal kicks (Giacobbo & Mapelli Reference Giacobbo and Mapelli2018, Reference Giacobbo and Mapelli2019a) and ComBinE use their own natal kick prescription (Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018). COMPAS and ComBinE use a distinct prescription for USSNe; StarTrack initially did not separately account for them, but Chruslinska et al. (Reference Chruslinska, Belczynski, Klencki and Benacquista2018) followed the model from Bray & Eldridge (Reference Bray and Eldridge2016), in which the natal kick is proportional to the ejecta mass. Similarly, MOBSE evolved from not explicitly accounting for USSNe (Giacobbo & Mapelli Reference Giacobbo and Mapelli2018) to treating USSNe as low-kick CCSNe (Giacobbo & Mapelli Reference Giacobbo and Mapelli2019a). Giacobbo & Mapelli (Reference Giacobbo and Mapelli2020) estimate the natal kick using the ejecta mass, similarly to Bray & Eldridge (Reference Bray and Eldridge2016).
In the context of mass transfer and CEEs, both crucial phases in DNS formation, the main differences appear in the definition of dynamical instability and in the treatment of stable mass transfer and CEEs.
For the definition of stability, MOBSE follows BSE (Hurley et al. Reference Hurley, Tout and Pols2002). StarTrack uses adiabatic mass loss indexes and a temperature threshold to identify stars with a fully convective envelope. COMPAS generally follows BSE and StarTrack but with significant amount of changes (see Vigna-Gómez et al. Reference Vigna-Gómez2018 and Appendix A). ComBinE use their own criteria based on a critical mass ratio (Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018).
COMPAS and MOBSE both follow BSE to solve for mass transfer episodes. In the default COMPAS models, non-accreted mass leaves the binary with the specific angular momentum of the accretor. StarTrack sets a fixed fraction of mass transfer as conservative, while any non-accreted mass leaves the system with the specific angular momentum of the orbit. Kruckow et al. (Reference Kruckow, Tauris, Langer, Kramer and Izzard2018) follows the approach from Soberman et al. (Reference Soberman, Phinney and van den Heuvel1997).
For the CEE, COMPAS closely follows StarTrack as implemented in Dominik et al. (Reference Dominik2012), particularly in the estimation of the binding energy structure parameter 
$\lambda$
(see Section 2.2). MOBSE follows Dewi & Tauris (Reference Dewi and Tauris2000) and Claeys et al. (Reference Claeys, Pols, Izzard, Vink and Verbunt2014) to determine this 
$\lambda$
parameter. Additionally, MOBSE frequently has highly efficient envelope ejection by allowing 
$\alpha > 1$
. Following Kruckow et al. (Reference Kruckow, Tauris, Langer, Sz´ecsi, Marchant and Podsiadlowski2016), ComBinE self-consistently uses their own stellar models to calculate the binding energy.
Most comparisons between population synthesis codes come from merger rates. However, the merger rate density is not the only prediction from population synthesis, and the full properties of observed populations will place stronger constraints on the physics (Barrett et al. Reference Barrett2018; Kruckow et al. Reference Kruckow, Tauris, Langer, Kramer and Izzard2018). The population of luminous red novae which will become DNSs do not seem to be an unequivocal additional constraint (see Section 4.2). Be X-ray binaries and short gamma-ray bursts must be included as additional constrains to the observed Galactic DNS population and the growing catalogue of DNS GW sources to fully dissect the origin and formation of DNSs.
5. Summary and conclusions
We carried out a rapid population synthesis study of a million massive binaries using COMPAS, finding 15 201 (unweighted) simulated systems which experience a CEE and eventually become a DNS. We present the key properties of the donor and binary star at the onset of the RLOF phase leading to the CEE. We provide an online catalogue of this synthesised population.
Some of our main results are
The CEEs that occur in DNS progenitors can be broadly divided into two types (description in Sections 3.1 and 4.1). Channel I, which accounts for 69% of all formed DNSs, involves a post-MS donor (the initially less massive star) with a NS companion. Channel II, which accounts for 14% of all formed DNSs, involves two giant stars, with mass ratio close to unity, in a double-core CE.
Close to 10% of the brightest luminous red nova transients, which have been previously associated with stellar mergers and CE ejections, are predicted to occur during binary evolution that leads to DNS formation. The progenitor of M101 OT2015-1 as reported in Blagorodnova et al. (Reference Blagorodnova2017) is somewhat similar to the pre-CEE properties of DNS-forming systems (see Section 4.2 and Figure 2).
We find that tidal circularisation timescales can be long compared to stellar radial growth timescales (see Figures 6 and 8), especially for rapidly evolving HG donors and/or donors with radiative envelopes experiencing only the less efficient dynamical tide rather than the more efficient equilibrium tide. This indicates that
$\sim 20\%$
–70% of binaries may not circularise prior to the onset of CEEs (see Sections 3.6 and 4.1). This finding suggests that the ensuing CE phases in these binaries may be distinct from those that have been previously considered. Future work is needed to determine the implication of these differences for the predicted formation rate and properties of DNSs.
One of the main goals of this study is to constrain the parameter space of interest for detailed evolutionary studies of CEEs. We hope that the results presented in this catalogue can inform choices of initial conditions for detailed hydrodynamical simulations and lead to an improved understanding of the complexities of dynamically unstable mass transfer and the subsequent CE phase. In particular, our present work highlights the roles of several uncertain processes that may be of crucial importance in DNS formation:
(i) Tidal dissipation in pre-CEE binary evolution;
(ii) Eccentric RLOF; and
(iii) The hydrodynamics of double-core CEEs.
Acknowledgements
We thank Poojan Agrawal, J.J. Eldridge, Mike Y.M. Lau, Orsola de Marco, Noam Soker, Simon Stevenson, and Eric Thrane for discussions. We thank Teresa Rebagliato for the graphic representation of the formation channels. A.V-G. acknowledges funding support from Consejo Nacional de Ciencia y Tecnología (CONACYT) and support by the Danish National Research Foundation (DNRF132). M.M. is grateful for support for this work provided by NASA through Einstein Postdoctoral Fellowship grant number PF6-170169 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. This material is based upon work supported by the National Science Foundation under Grant No. 1909203. S.dM. and S.J. acknowledge funding by the Netherlands Organization for Scientific Research (NWO) as part of the Vidi research programme BinWaves with project number 639.042.728 and the European Union’s Horizon 2020 research and innovation programme from the European Research Council (ERC, grant agreement No. 715063). I.M. is a recipient of the Australian Research Council Future Fellowship FT190100574.
A. Population synthesis: details of the COMPAS set-up
We present a list of the initial values and default settings used for this study in Table A1 in order to be able to emulate them with other population synthesis codes. References have been added where needed in order to justify our assumptions. Some of these assumptions are described in Section 2.
Initial values and default settings of the population synthesis model simulations with COMPAS, including pertinent references.
