Common envelope episodes that lead to double neutron star formation

Abstract

Close double neutron stars (DNSs) have been observed as Galactic radio pulsars, while their mergers have been detected as gamma-ray bursts and gravitational wave sources. They are believed to have experienced at least one common envelope episode (CEE) during their evolution prior to DNS formation. In the last decades, there have been numerous efforts to understand the details of the common envelope (CE) phase, but its computational modelling remains challenging. We present and discuss the properties of the donor and the binary at the onset of the Roche lobe overflow (RLOF) leading to these CEEs as predicted by rapid binary population synthesis models. These properties can be used as initial conditions for detailed simulations of the CE phase. There are three distinctive populations, classified by the evolutionary stage of the donor at the moment of the onset of the RLOF: giant donors with fully convective envelopes, cool donors with partially convective envelopes, and hot donors with radiative envelopes. We also estimate that, for standard assumptions, tides would not circularise a large fraction of these systems by the onset of RLOF. This makes the study and understanding of eccentric mass-transferring systems relevant for DNS populations.

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 1976). 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 1976). Most evolutionary pathways leading to close compact binaries are expected to involve at least one CEE (Ivanova et al. 2013a).

While CEEs are frequently invoked as a fundamental part of binary evolution, the detailed physics remain poorly understood (Iben & Livio 1993; Ivanova et al. 2013b; Paczynski 1976). There have been efforts in modelling and understanding the phase through hydrodynamic simulations, using Eulerian adaptive mesh refinement (Chamandy et al. 2018; De et al. 2019; Iaconi et al. 2017, 2018; Li et al. 2020; L´opez-C´amara et al. 2019; MacLeod & Ramirez-Ruiz 2015; MacLeod et al. 2017b; Passy et al. 2012; Ricker & Taam 2008, 2012; Sandquist et al. 1998; Shiber et al. 2019; Staff et al. 2016), moving meshes (Ohlmann et al. 2016, 2017; Prust & Chang 2019), smoothed particle (Ivanova & Nandez 2016; Lombardi et al. 2006; Nandez, Ivanova, & Lombardi 2015; Nandez & Ivanova 2016; Passy et al. 2012; Rasio & Livio 1996; Reichardt et al. 2019), particle-in-cell (Livio & Soker 1988), and general relativistic (Cruz-Osorio & Rezzolla 2020) methods. Other approaches pursue detailed stellar modelling (Clayton et al. 2017; Dewi & Tauris 2000; Fragos et al. 2019; Klencki et al. 2020; Kruckow et al. 2016) or binary population synthesis (e.g., Andrews et al. 2015; Dewi, Podsiadlowski, & Sena 2006; Kruckow et al. 2018; Nelemans et al. 2000; Tauris & Bailes 1996; Vigna-Gómez et al. 2018). 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. 2018; Chruslinska et al. 2017, 2018; Dominik et al. 2012), MOBSE (Giacobbo & Mapelli 2018, 2019a,b, 2020), ComBinE (Kruckow et al. 2018), and Compact Object Mergers: Population Astrophysics and Statistics (COMPAS) (Chattopadhyay et al. 2020; Vigna-Gómez et al. 2018) 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. 2015; Belczynski, Kalogera, & Bulik 2002; Belczynski et al. 2018; Bhattacharya & van den Heuvel 1991; Dewi & Pols 2003; Dewi, Podsiadlowski, & Pols 2005; Ivanova et al. 2003; Kruckow et al. 2018; Tauris & van den Heuvel 2006; Tauris et al. 2017; Vigna-Gómez et al. 2018). Dominik et al. (2012) 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 cataloguea.

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. 2017; Howitt et al. 2020; Ivanova et al. 2013b; Pastorello et al. 2019; MacLeod et al. 2017a), and Be X-ray binaries (Vinciguerra et al. 2020).

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 suiteb (Barrett et al. 2018; Neijssel et al. 2019; Stevenson et al. 2017; Vigna-Gómez et al. 2018). 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. 2009). 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 2001) 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. 2012). 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 1924). 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. 2019). 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. (2018). 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. 1. We approximate the Roche lobe radius following the fitting formula provided by Eggleton (1983) in the form:

   (1)\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}
   where $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. 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)\begin{equation} \zeta_{\textrm{i}}=\frac{d\log R_{\textrm{i}}}{d\log m_{\textrm{donor}}}, \end{equation}
   where the subscript ‘i’ represents either the mass-radius exponent for the donor ( $\zeta_{\textrm{donor}}$ ) or for the Roche lobe ( $\zeta_{\textrm{RL}}$ ). We assume that
   (3)\begin{equation} \zeta_{\textrm{donor}} < \zeta_{\textrm{RL}} \end{equation}
   leads to a CEE. Inspired by Ge et al. (2015), for MS donors, we assume $\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. (1997). All mass transfer episodes from stripped post-helium ignition stars, that is, case BB mass transfer (Delgado & Thomas 1981; Dewi et al. 2002; Dewi & Pols 2003) neutron star (NS) are assumed to be dynamically stable. For more details and discussion, see Tauris et al. (2015) and Vigna-Gómez et al. (2018).

3. 3. We deviate from Stevenson et al. (2017) and Vigna-Gómez et al. (2018) 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. 4. We follow Kool (1990) in the parameterisation of the binding energy ( $E_{\textrm{bind}}$ ) of the donor star’s envelope ( $m_{\textrm{donor,env}}$ ) given as:

   (4)\begin{equation} E_{\textrm{bind}}=\frac{-Gm_{\textrm{donor}}m_{\textrm{donor,env}}}{\lambda R_{\textrm{donor}}}, \end{equation}
   where G is the gravitational constant and $\lambda$ is a numerical factor that parameterises the binding energy.

5. 5. For the value of the $\lambda$ parameter, we follow the fitting formulae from detailed stellar models as calculated by Xu & Li (2010a,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. 2012). 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. 6. We use the $\alpha\lambda$ -formalism (de Kool 1990; Webbink 1984) to determine the post-CEE orbit, with $\alpha=1$ in all of our CEE.

7. 7. We use the Fryer et al. (2012) delayed supernova (SN) remnant mass prescription, which was the preferred model from Vigna-Gómez et al. (2018). 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. 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 1977). Meanwhile, the dynamical tide (Zahn 1975) 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 1998; Moe & Kratter 2018). 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 2008).

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. 2016 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 1999a,b). The high-eccentricity regime, previously studied in the parabolic (and chaotic) limit (Mardling 1995a,b), has recently being revisited in the context of both dynamical (Vick & Lai 2018) and equilibrium tides (Vick & Lai 2020). 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 2008).

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 (1981) and implemented by Hurley, Tout, & Pols (2002), although this may not be accurate for high-eccentricity systems (but see Vick & Lai 2020). 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 (1977), as used in Hurley et al. (2002).

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

(5)\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

(6)\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 (1981). The structure of the tidally deformed star is parameterised by k, which is the apsidal motion constant (Lecar, Wheeler, & McKee 1976) and the intrinsic tidal timescale ( $\tau_{\textrm{tide}}$ ), usually associated with viscous dissipation (Zahn 1977). We follow Hurley et al. (2002) 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

(7)\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 (1977), we can write the synchronisation and circularisation timescales for the dynamical tide as

(8)\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

(9)\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. (2002) from the values given by Zahn (1975), 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 (2000) to the detailed stellar models from Pols et al. (1998).

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 2008) and implementation (Hurley et al. 2002). Siess et al. (2013) 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. 2017). The calculation of k and $\tau_{\textrm{tide}}$ follows Hurley et al. (2000) 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 2019).

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. 2012; Giacobbo & Mapelli 2018; Neijssel et al. 2019; Vigna-Gómez et al. 2018).

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. (2018). 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.

Figure 1. 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. 2017 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 1991; Tauris & van den Heuvel 2006; Tauris et al. 2017). Following Vigna-Gómez et al. (2018), 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. 1. A post-MS primary engages in stable mass transfer onto a MS secondary.

2. 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. 3. The secondary evolves off the MS, expanding and engaging in a CEE with the NS accretor.

4. 4. After successfully ejecting the envelope, and hardening the orbit, the secondary becomes a naked helium star.

5. 5. The stripped post-helium-burning secondary engages in highly non-conservative stable (case BB) mass transfer onto the NS companion.

6. 6. After being stripped of its helium envelope, the ultra-stripped secondary (Tauris et al. 2013, 2015) 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. 2012). 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. 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 1995; Dewi et al. 2006; Justham, Podsiadlowski, & Han 2011). 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. (2006). During this double-core CEE, both stars are stripped and become naked helium stars.

2. 2. The stripped post-helium-burning primary engages in stable (case BB) mass transfer onto a stripped helium burning secondary.

3. 3. The primary, now a naked metal star, explodes in a SN and becomes a NS.

4. 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. 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 2001), comprise less than 2% of the DNS population.

3.2. Comparison with Vigna-Gómez et al. (2018)

This work generally uses similar assumptions and physics parameterisations as the preferred model of Vigna-Gómez et al. (2018), including the Fryer et al. (2012) 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. (2018). 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. (2018). Here, we estimate the formation rate of all (merging) DNS to be 85 (60) $\rm Gpc^{-3}/yr$ . Vigna-Gómez et al. (2018). 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. (2000)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.

Figure 2. 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. 2017) with a star symbol. The solid black lines indicate ZAMS and TAMS loci for a grid of SSE models (Hurley et al. 2000) 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. (2016) and the dashed grey line shows the stellar track from Pols et al. (1998, 2009). The dash-dotted light blue and solid green lines show how fitting formulae from Hurley et al. (2000) 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. 2000). 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. (2000) 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. (2016) also display a bifurcation point, but models and interpolation are smoother than those in Hurley et al. (2000).

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. (2000). Donors which engage in a CEE leading to DNS formation can be in the HG (4%), GB (7%), CHeB (59%), or EAGB (30%) phase.

Table 1. 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.

Figure 3. 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.

Figure 4. 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 1977 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}$ .

Figure 5. 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 (1991); Tauris & van den Heuvel (2006); Tauris et al. (2017) and references therein. Channel II involves a double-core CE between two post-MS stars. A similar channel has been proposed by Brown (1995) and Dewi et al. (2006), 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. (2006). 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).

Table 2. Distinct DNS subpopulations as described in Section 4.1 and presented in Figure 7.

Figure 6. 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.

Figure 7. 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. (2008) 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. (2020) 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.

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. (2002) 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. (2008) 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.

Figure 9. 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. (2000) do not contain explicit information about the moment when this transition occurs. Hurley et al. (2002) 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. (2008) 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. (2000). 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. (2017). This event is similar to other luminous red novae associated with CEEs (Ivanova et al. 2013b). Following the discovery of M101 OT2015-1, archival photometric data from earlier epochs were found. Blagorodnova et al. (2017) 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).

Figure 10. 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. (2017) 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. (2020) 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 1996; Portegies Zwart & Yungelson 1998; Toonen, Nelemans, & Portegies Zwart 2012), STARTRACK (Belczynski et al. 2002, 2008), binary stellar evolution (BSE) (Hurley et al. 2002), the Brussels code (De Donder & Vanbeveren 2004), COMPAS (Stevenson et al. 2017; Vigna-Gómez et al. 2018), ComBinE (Kruckow et al. 2018), and customised software based on them all assume that RLOF commences in circular binaries. Detailed simulations, such as those of Passy et al. (2012), MacLeod, Ostriker, & Stone (2018), and others, often make the assumption of an initially circular orbit (but see Staff et al. 2016, 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 1983, 1984; Sepinsky et al. 2007, 2009, 2010; Dosopoulou & Kalogera 2016a,b).

The analysis of Sepinsky et al. (2007) 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 (2016b) 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 (2019) noted that evolution towards circularisation from Dosopoulou & Kalogera (2016b) 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 (2016b). 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. 2015). Mass transfer onto a NS companion is highly non-conservative, almost in the fully non-conservative limit (Tauris et al. 2015).

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 2017; Church et al. 2009; Lajoie & Sills 2011; Regös, Bailey, & Mardling 2005; Staff et al. 2016; van der Helm, Portegies Zwart, & Pols 2016).

Staff et al. (2016) 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 (2019) 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 2018, see also Schrøder et al. 2020). The former results in a less luminous transient, comparable in energetic to luminous red novae (Kashi & Soker 2016). The interaction is driven by jets which might enhance mass loss at periastron passages, keeping the system eccentric (Kashi & Soker 2018). These jets might, in some cases, prevent CEEs (Shiber & Soker 2018).

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 1999). Eccentric ( $e\lesssim 0.2$ ) MS White Dwarf binaries which are believed to have experienced RLOF are not rare (Kawahara et al. 2018; Masuda et al. 2019; Vos et al. 2013). Jayasinghe et al. (2019) 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. (2019) will engage in RLOF at some later point.

Sirius (Gatewood & Gatewood 1978; van den Bos 1960) 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 (2008) 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. (2008), Saladino & Pols (2019) 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. (2007) 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 (2009) discusses $\gamma^2$ Velorum as post-mass transfer binary system. In that work, Eldridge (2009) 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. 2014; Hurley et al. 2002).

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. 2010; Hjellming & Webbink 1987; Ivanova & Taam 2004). 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 2002). 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. (2017) 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. (2019) 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. (1996), Hurley et al. (2002), and references therein.

Witte & Savonije (1999a,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 1990; Levato et al. 1987; Kobulnicky et al. 2014; Sana et al. 2012 as presented in Figure 3 of Moe & Di Stefano 2017). 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. (2018) 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. 2018 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. 2012; Kruckow et al. 2018; Neijssel et al. 2019). 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 2014; Podsiadlowski, Joss, & Hsu 1992; Vigna-Gómez et al. 2019).

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 2019 but see, e.g., Ye et al. 2020). Meanwhile, Kozai–Lidov oscillations (Kozai 1962; Lidov 1962) in hierarchical triple systems can drive up the eccentricity of the inner binary (see Naoz et al. 2016 for a review) and contribute to the formation of merging DNSs (Hamers & Thompson 2019). 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. 2019) based on a single detection of GW170817 with a flat-in-rate prior (Abbott et al. 2017). 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. 2020).

Vigna-Gómez et al. (2018) predicted a DNS merger rate of ${\approx}280$ Gpc⁻³/yr. After several changes in COMPAS (see Section 3.2), the DNS merger rate is now ${\approx}60$ Gpc⁻³/yr. This is more in line with the reported rates from other rapid population synthesis codes (Giacobbo, Mapelli, & Spera 2018; Giacobbo & Mapelli 2019b; Kruckow et al. 2018). Chruslinska et al. (2018) 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 (2020) proposed that a revised natal kick prescription might prove important to reconcile DCO merger rates, as suggested by Chruslinska et al. (2018) (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 (2015). They find that the rate can increase (decrease) by a factor of a few for plausible shallower (steeper) initial mass functions. Schneider et al. (2018a, see also Farr & Mandel 2018; Schneider et al. 2018b) 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 (1955) value. While this level of fluctuation in the initial mass function may prove insufficient to resolve the DNS rate discrepancy (Belczynski et al. 2018), 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 (2020). This can be converted to a volumetric rate in units of Gpc⁻³/yr by multiplying the rate per Myr per Milky Way equivalent Galaxy by a factor of $\sim 10$ (Abadie et al. 2010). 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. 2010), and all rate intervals are broad.

Lau et al. (2020) use a synthesised DNS population to predict that a 4-yr Laser Interferometer Space Antenna (LISA) mission (Amaro-Seoane et al. 2017; Baker et al. 2019) will detect 35 Galactic DNSs. Andrews et al. (2020) 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 (2015) presents 69 systems. From those 69, 44 have an observed orbital period ( $P_{\textrm{orb}}<520$ d). Vinciguerra et al. (2020) 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. (2020) 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. (2012) for the remnant mass prescription of compact objects, using either the rapid or delayed variation, while ComBinE use their own (Kruckow et al. 2018). The natal kick distribution varies not only between codes but also between papers. COMPAS initially followed StarTrack (e.g., Belczynski et al. 2008, 2018) with a no-kick model for electron-capture SN (ECSN) (Stevenson et al. 2017) 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 2018, 2019a) and ComBinE use their own natal kick prescription (Kruckow et al. 2018). COMPAS and ComBinE use a distinct prescription for USSNe; StarTrack initially did not separately account for them, but Chruslinska et al. (2018) followed the model from Bray & Eldridge (2016), in which the natal kick is proportional to the ejecta mass. Similarly, MOBSE evolved from not explicitly accounting for USSNe (Giacobbo & Mapelli 2018) to treating USSNe as low-kick CCSNe (Giacobbo & Mapelli 2019a). Giacobbo & Mapelli (2020) estimate the natal kick using the ejecta mass, similarly to Bray & Eldridge (2016).

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. 2002). 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. 2018 and Appendix A). ComBinE use their own criteria based on a critical mass ratio (Kruckow et al. 2018).

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. (2018) follows the approach from Soberman et al. (1997).

For the CEE, COMPAS closely follows StarTrack as implemented in Dominik et al. (2012), particularly in the estimation of the binding energy structure parameter $\lambda$ (see Section 2.2). MOBSE follows Dewi & Tauris (2000) and Claeys et al. (2014) to determine this $\lambda$ parameter. Additionally, MOBSE frequently has highly efficient envelope ejection by allowing $\alpha > 1$ . Following Kruckow et al. (2016), 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. 2018; Kruckow et al. 2018). 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. (2017) 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:

1. (i) Tidal dissipation in pre-CEE binary evolution;

2. (ii) Eccentric RLOF; and

3. (iii) The hydrodynamics of double-core CEEs.