© The Authors 2025

1. Introduction

The discovery of gravitational waves (GW) from merging binary compact objects (Abbott et al. 2016) has revolutionised observational astronomy. The observations made by the Advanced LIGO (LIGO Scientific Collaboration 2015) and Advanced Virgo (Acernese et al. 2015) detectors have led to major breakthroughs in the study of compact objects and their stellar progenitors and opened new avenues for fundamental physics and cosmology (Abbott et al. 2023a). While the vast majority of the detected sources are binary black holes (BBHs), two binary neutron star (BNS) mergers have been observed. The first detection, GW 170817 (Abbott et al. 2017a,b), was a multi-messenger source, having been observed in the optical (Coulter et al. 2017), infrared, X-ray, gamma-ray (Goldstein et al. 2017; Savchenko et al. 2017), and radio (Mooley et al. 2018; Ghirlanda et al. 2019). These observations confirmed the long-standing hypothesis that some short gamma-ray bursts originate from a jet launched during the BNS merger and allowed for the study of the geometry of the jet. In addition, observations in infrared, optical, and ultraviolet confirmed the presence of a kilonova. This phenomenon consists of the thermal radiation from the radioactive decay of massive elements produced in the optically thick BNS merger ejecta by rapid neutron capture (r-process). This kilonova observation confirmed that BNS mergers are sites of r-process nucleosynthesis (see e.g. Villar et al. 2017). The second BNS detection (GW 190425), while lacking observed electromagnetic counterparts, is also remarkable in that the total mass of the system (3.4 M_⊙) is significantly larger (Abbott et al. 2020) than any of the other known BNSs (Galactic pulsars as well as GW 170817). These two detections also allowed the merger rate in the local Universe to be estimated, albeit with a large uncertainty, as 10 − 1700 Gpc⁻³ ⋅ yr⁻¹ (Abbott et al. 2023b). It is interesting to note that Galactic binary pulsars and BNS systems observed via their GW emission probe different evolutionary stages of the same underlying population, ranging from large separation for binary pulsars to small separation (due to GW emission) and merger for BNSs detected with GWs.

While the number of known BNS systems is currently relatively small, it is expected to grow in the very near future. Based on current population models, at design sensitivity the LIGO-Virgo-Kagra network is expected to observe up to about 100 systems per year (Kiendrebeogo et al. 2023). Moreover, the third-generation GW detectors (Einstein Telescope and Cosmic Explorer) planned for the 2030s are expected to detect up to 10⁵ BNS mergers per year (Maggiore et al. 2020; Evans et al. 2021; Branchesi et al. 2023).

The formation scenarios of binary compact objects, and in particular BNSs, have been extensively studied using methods with different degrees of approximations. In particular, population synthesis codes (e.g. Portegies Zwart & Verbunt 1996; Hurley et al. 2002; Belczynski et al. 2002; Spera et al. 2015; Eldridge & Stanway 2016; Stevenson et al. 2017; Giacobbo et al. 2018; Kruckow et al. 2018; Breivik et al. 2020; Riley et al. 2022; Fragos et al. 2023; Iorio et al. 2023) have become a standard tool in compact binary studies thanks to their computational efficiency, and they have been used to estimate the BNS merger rates (see references in Mandel & Broekgaarden 2022). Most of these codes follow stellar evolution using detailed simulations for single stars, assuming that they are not affected by the presence of the companion. In contrast, BPASS (Eldridge & Stanway 2016) and POSYDON (Fragos et al. 2023) rely on detailed binary evolution simulations. Stellar properties at each time step are then computed using either fitting formulae (e.g. COSMIC, Breivik et al. 2020; COMPAS, Riley et al. 2022; StarTrack, Belczynski et al. 2002) or interpolation on pre-computed grids (e.g. SEVN; Iorio et al. 2023; POSYDON). In population synthesis codes based on single-star evolution, the processes related to binarity (e.g. mass exchanges, orbital evolution) are computed using phenomenological prescriptions. As most BNSs are likely to be formed by isolated binaries, this method is particularly useful since it allows the cosmological populations of compact binaries to be modelled and large portions of the parameter space related to binary evolution to be probed. In addition, several codes have been developed for the study of dynamical effects in dense stellar systems (e.g. Giersz et al. 2013; Kremer et al. 2020). While the details of different population synthesis codes vary, several general conclusions on BNS populations have emerged.

First of all, many recent studies have found that most BNSs are expected to form from isolated binary massive stars, rather than through dynamical interactions in dense environments (Belczynski et al. 2018; Zevin et al. 2019; Ye et al. 2020; Fragione & Banerjee 2020). Indeed, while three- and four-body interactions can drive BBHs to merge in dense stellar clusters, this is not the case for BNSs. The reason for this difference is that BBHs, being more massive, dominate the cores of stellar clusters, thereby preventing mass segregation of neutron stars (NSs). This effect greatly reduces the number of dynamical interactions NSs can undergo. Additionally, because NSs are less massive, they are more likely to be ejected out of the clusters upon formation due to their natal kicks.

The second conclusion is that within the isolated formation channel, a specific evolutionary track is favoured. It features a binary massive star at zero age main sequence (ZAMS), a phase of common envelope (CE) between a giant and a NS, and a later phase of mass transfer between the stripped giant star and the NS (case BB mass transfer, e.g. Tauris et al. 2015). Both mass transfer episodes shrink the orbit considerably so that the resulting binary compact object is more likely to merge within the Hubble time t_Hubble = 13.7 Gyr. Moreover, the secondary explodes as an ultra-stripped supernova (SN), which leads to a smaller natal kick compared to a standard core-collapse SN, thereby increasing the probability that the system remains gravitationally bound after the second SN.

More generally, several studies have pointed out that the natal kicks of core-collapse SNe that form NSs should be weaker than those that form BHs, of the order of 20 km ⋅ s⁻¹ (e.g. Giacobbo & Mapelli 2018). Conversely, BH natal kicks are typically one order of magnitude larger (Hobbs et al. 2005), though some more recent models scale them down by a factor 1.4 M_⊙/M_BH so that NSs and BHs receive the same momentum after the collapse (e.g. Riley et al. 2022). Fallback accretion may further reduce the BH natal kick (Fryer et al. 2012). Low NS kicks are naturally expected in the case of ultra-stripped SNe, electron capture SNe, and accretion-induced collapse (AIC) of a white dwarf (WD) into a NS. Such low kicks in AIC events may increase the survival rate of BNS progenitors (Chruslinska et al. 2018). As a result, AIC could play an important role in the formation of BNSs. We briefly summarise this process in the following paragraph.

Accretion-induced collapse occurs when a heavy (O/Ne/Mg or Si/O) white dwarf accretes matter from its stellar companion with a slow enough rate (≲ 10⁻⁹ M_⊙ · yr⁻¹), leading the WD to collapse and form a NS. In general, AIC can take place either in a double-degenerate case (WD-WD merger) or a single-degenerate case (WD accreting from a companion star that overflows its Roche lobe). Ruiter et al. (2019) suggest that double-degenerate AIC events operating in WD binary systems could be an important step in forming merging BNSs observed with ground-based GW interferometers. In this work, we show the possibility that single-degenerate AIC events can also contribute to the population of merging BNSs at high redshifts. AIC events are expected to produce faint optical transients, but these would be short-lived (lasting from a few days to a week) and under-luminous with respect to regular SNIa (Longo Micchi et al. 2023). Such events are therefore difficult to observe, and there have been no detections reported to date, although some indirect evidence exists (see e.g. Wang & Liu 2020, for a review).

Previous studies of BNS populations have also explored the time delays Δt_delay between the formation of the stellar binary (ZAMS) and the merger of the BNS. This delay is defined as the sum of the stellar evolution part (the time for the massive stellar binary to form a compact object binary), Δt_evol, and the orbital evolution due to the emission of GWs (orbital decay of a compact object up to coalescence), Δt_GW. Usually, a single distribution of the time delays is assumed for an entire cosmological population, and it is found that the distribution closely follows a power law, P(Δt_delay) ∝ Δt_delay⁻¹ (e.g. Chruslinska et al. 2018; Neijssel et al. 2019). In fact, this is the expected distribution if the delays are dominated by the GW phase (i.e. assuming negligible stellar evolution time) and if the distribution of the initial semi-major axes a_i is log-uniform since the GW delay scales as Δt_GW ∝ a_i⁴. However, the question arises as to whether this distribution is universal, that is, whether it applies to all the evolutionary tracks of BNS progenitors at all epochs or if it depends on the stellar evolution parameters.

In this work, we explore the evolutionary tracks of BNS progenitors using the population synthesis code COSMIC (Breivik et al. 2020) and study in detail the distribution of delay times through cosmic history. Throughout the paper, we focus on systems that merge within the age of the Universe. Our population synthesis parameter choices, metallicity sampling protocol, and procedure to identify the evolutionary tracks and the merger rate density calculation are explained in Sect. 2. With our assumptions on stellar evolution parameters, we identified three evolutionary tracks: a standard one involving case BB mass transfer; a track involving two stars with almost identical ZAMS masses, which leads to synchronised evolution and a phase of CE with both stars in their giant phase; and a track that involves the AIC of one of the progenitors. These tracks are discussed in Sect. 3.1. We then focus on the dependence on metallicity and explore in detail the properties of the progenitor stellar binaries that contribute to each track in Sect. 3.2. We discuss the properties of the BNS populations in each track in Sect. 3.3. We find that the time delay distributions at a given metallicity do not follow a simple power law but instead present a complex structure that depends on the evolutionary track. Having identified the dominant evolutionary tracks at each metallicity, we explore the properties of the resulting cosmological population in Sect. 3.4. In particular, we find that the AIC evolutionary track has very short delay times Δt_GW and is dominant at high redshift. However at low redshift, the equal-mass evolutionary track is dominant, and it is the one expected to contribute to the population observed with GW detectors. We discuss some interesting applications and consequences of our results as well as some caveats to our analysis in Sect. 4.

2. Modelling binary neutron star populations

2.1. Population synthesis code

In this work, we use the COSMIC population synthesis code (Breivik et al. 2020) to simulate a population of BNSs. COSMIC implements stellar evolution using Hurley et al. (2000) and binary interactions using the binary stellar evolution (BSE) code from Hurley et al. (2002), with several modifications to account for recent updates to binary evolution. We use the default parameters of version 3.4.0, with a few exceptions that are discussed in this section. The full parameter list used in this study can be found in Appendix A. We note that the goal of this study is not to perform a systematic sampling of the possible parameter values but rather to focus on one physically motivated model and focus on the evolutionary tracks that lead to the formation of the BNSs that we describe in Sect. 3.1.

We assumed standard properties of the binaries at ZAMS. The masses of the primary follow a broken power-law distribution (Kroupa 2001), and mass ratios are uniformly distributed (e.g. Sana et al. 2012). We followed the distributions of initial semi-major axes and eccentricities introduced by Sana et al. (2012).

Stellar winds are treated using recent corrections to BSE discussed in Breivik et al. (2020) (corresponding to windflag = 3). In this study, we do not introduce Eddington-limited winds as proposed by Gräfener & Hamann (2008), where stronger winds are less affected by metallicity (see also Giacobbo et al. 2018). The average accretion rate per orbit onto the companion is assumed to follow the Bondi & Hoyle (1944) mechanism and is therefore limited. In particular, most mass transfer events are non-conservative. Thus, any event of AIC is only possible following a phase of stable mass transfer onto a WD companion. Stellar winds mostly impact binary evolution because of the single-star loss of mass. Wind velocities are assumed to depend on the stellar type (Belczynski et al. 2008).

Due to the envelope expansion in the late stages of stellar evolution, the stellar radius can exceed its Roche-lobe radius (R_L), which is always computed at periastron and assuming the expression from Eggleton (1983). When mass transfer starts, the orbit is assumed to instantaneously circularise at periastron. Mass transfer can be stable (Roche lobe overflow) or unstable, in which case this leads to a CE (see e.g. the review by Ivanova et al. 2013). The treatment of the boundary between stable and unstable mass transfer rates depends on the stellar type of (only) the donor star and on the mass ratio q = M_don/M_comp. If q > q_crit, mass transfer is dynamically unstable and a phase of CE follows. We use the values for q_crit proposed by Neijssel et al. (2019) and also used in e.g. COMPAS (qcflag = 5). In particular, with these values, mass transfers from stripped donors are always dynamically stable.

In the isolated binary formation scenario, the CE phase allows the orbital separation to shrink enough for the compact object binary formed after the stellar explosions to merge within t_Hubble. COSMIC treats CE in a simple parametric way using the alpha lambda formalism (Webbink 1984), where α_CE describes the fraction of the orbital energy that is transferred to the envelope, leading to its expansion and possible ejection; and λ (of order unity) describes the stiffness of the envelope density profile. In this study, we fix α_CE = 1 (we do not include external energy terms as discussed in Fragos et al. 2019; Santoliquido et al. 2021); and a variable λ that depends on the stellar type (see Claeys et al. 2014). We assume that stellar companions without a clear core-envelope boundary automatically lead to a merger during the CE phase (cemergeflag = 1), as introduced by Belczynski et al. (2008). We do not deduce the final core mass of the donor star in the case of case BB mass transfer using the expressions for the rates from Tauris et al. (2015), and instead assume the expression from Hurley et al. (2002) (cehestarflag = 0).

Natal kicks are given to the compact objects that form following SNe, mostly because of asymmetries in the ejected material during the explosion (see e.g. Janka & Mueller 1994; Wongwathanarat et al. 2015). For core-collapse SNe, we assume that the natal kicks follow a Maxwellian distribution with a kick velocity dispersion σ_k = 265 km ⋅ s⁻¹ (Hobbs et al. 2005). For electron-capture SNe, ultra-stripped SNe, and AIC events, we assumed that the kicks are reduced and have a dispersion σ_k, low = 20 km ⋅ s⁻¹. This parameter is especially important in this study, as most evolutionary tracks involve one or several of these events with lower kicks and therefore higher probabilities for the binary to survive the explosions. We do not assume that the kick velocity is affected by the ejected mass or the remnant mass as, for example, Giacobbo & Mapelli (2020) (kickflag = 0). We infer the remnant mass following the rapid mechanism for the SN explosion (Fryer et al. 2012), with updates from Giacobbo & Mapelli (2020) (remnantflag = 3). This leads to a mass gap between NSs and BHs, contrary to the delayed mechanism.

Finally, the Solar metallicity is set to Z_⊙ = 0.014 (Asplund et al. 2009).

2.2. Sampling protocol

We create a grid of metallicities ranging from Z = 9.5 × 10⁻⁵ (6.8 × 10⁻³ Z_⊙, the minimum metallicity allowed in BSE (Hurley et al. 2002) and thus in COSMIC) to Z = 0.028 (2 Z_⊙). In practice, this last metallicity bin is not used in our calculation of the merger rate, and the details of stellar evolution at super-solar metallicities are more uncertain. The metallicity grid we use is thus 0.000095; 0.00014; 0.00021; 0.0003; 0.00044; 0.00065; 0.00095; 0.0014; 0.0021; 0.003; 0.0044; 0.0065; 0.0095; 0.014(; 0.028), i.e. 14 (15) metallicity values.

For each of these metallicities, we sample and evolve 9.55 × 10⁹ binaries. We do not use the match feature of COSMIC that automatically stops the sampling once the properties of the BNS masses and/or semi-major axes and/or eccentricities have converged to a stable distribution (for a more complete description, see Breivik et al. 2020). In our case, we decide to use a single random seed across all metallicities, meaning that the initial population is the same for all metallicities. For a given binary, stochastic parameters, such as SN kick intensity and orientations, are also seeded, which means that a direct comparison between individual binary systems across the metallicity bins is possible. A potential downside of our approach is that in some metallicity bins, the total number of BNS systems may be quite low (∼10⁴), and thus the studied distributions are slightly biased by the random seed. For a more accurate simulation, the random seed should also be varied across metallicities, but this is not expected to impact the results provided the BNS sample is large enough.

We focus on all binaries that produce BNS systems that remain gravitationally bound. Some of them have initial orbital properties (semi-major axis and eccentricity) that prevent them from merging within the Hubble time. We mark them as non-merging systems and remove them from our study.

2.3. Classification of binary evolutionary tracks

In order to define evolutionary tracks, we design an automated classification scheme based on the output of COSMIC. We note that analogous procedures were employed by several authors in the context of BBH and NSBH systems (see e.g. Vigna-Gómez et al. 2018; Neijssel et al. 2019; Bavera et al. 2021; Broekgaarden et al. 2021, 2022; Iorio et al. 2023). In our case, for each simulation at a given metallicity, we store the information on the binaries that produce gravitationally bound BNS systems, including all the stellar properties and binary properties at key stages of stellar evolution.

The first stage is to identify evolutionary sequences, which we define as the sequence of different stellar types for both the primary and the secondary (see Table 1 for the list of stellar types used in COSMIC). An illustration of this procedure is presented in Fig. 1, calculated at solar metallicity. The vertical (primary stellar type) and horizontal (secondary stellar type) axes refer to evolutionary sequences of the primary and secondary, respectively. Each line (column) for the primary (secondary) corresponds to a particular sequence of stellar types as the evolution of the binary proceeds. For example, the sequence 1- 1 2 3 4 5 8 9 13 corresponds to primary stars which evolve off the main sequence (1-2-3) until they start burning helium in their core (4). After a brief time on the asymptotic giant branch (5), their helium core is stripped of its envelope and the naked He star evolves (8-9) until its explosion, forming a NS (13). All primary evolutionary sequences containing 12 correspond to systems forming through AIC. The intersections (in colour) mark the combinations that allow for a formation of at least a BNS that merges within the age of the Universe. The colour coding shows the number of systems in each evolutionary sequence. While the total number of sequences that leads to a BNS is quite large, only a few of them are sufficiently frequent. In the example shown here, the 12 most frequent evolutionary sequences make up to 95% of the population of merging BNSs.

Fig. 1. Illustration of our classification scheme. Shown is the number of BNS systems for each sequence of stellar types of the primary (rows) and of the secondary (columns) at Z = Z⊙ = 0.014. The combination of primary and secondary evolutionary sequences creates 46 evolutionary sequences, but 95% of the systems are found in 12 of them. The sequences of numbers on the y axis and on the x axis are the series of stellar types (see Table 1) for the primary and the secondary, respectively. Only systems with Δtdelay < tHubble are shown.

This method does not guarantee a priori that systems with the same evolutionary sequence have similar binary evolution history, since we did not account for binary interactions (such as CE), nor the order of specific events (such as SNe) in this classification. However, we verified that for all the evolutionary sequences defined with this method, the details of binary stellar evolution for all the binaries are indeed very similar.

The dominant evolutionary sequences are not necessarily the same at all metallicities. We therefore repeat the same process over the populations at different metallicities to find all the dominant sequences. Overall, we extracted the 13 evolutionary sequences that make up most of the population of BNSs across the different metallicity bins and analysed each of them individually. We found that these 13 evolutionary tracks can in fact be grouped in only three categories with different binary histories, that we describe in more detail in the Sects. 3.1.1–3.1.3 below. In practice, we therefore leave a small fraction (∼5%) of systems unclassified, but this doesn’t affect the results presented here.

2.4. Merger rate density

We use the procedure described in Lehoucq et al. (2023) to combine the populations from all the metallicity bins and simulate a realistic astrophysical population. Here we briefly summarise the main ingredients of the calculation.

Redshift and time are related by the relation

$$\begin{array}{c}\hfill \frac{\mathit{dt}}{\mathit{dz}}=\frac{H_0^{-1}}{(1+z)\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}},\end{array}$$(1)

where H₀ = 68 km ⋅ s⁻¹ ⋅ Mpc⁻¹ is the Hubble constant, Ω_Λ = 0.69 is the fraction of dark energy in the total Universe energy density and Ω_m = 0.31 is the fraction of energy in matter (dark matter and baryonic matter). We also define the useful quantities Ω_b = 0.045, the baryonic fraction of energy density and h₀ = 0.68, the reduced Hubble constant (dimensionless). We assume here the cosmological values of Planck Collaboration VI (2020).

We used the following functional form for the star formation rate (SFR) from Springel & Hernquist (2003):

$$\begin{array}{c}\hfill \psi (z)=\nu \frac{aexp[b(z-z_{\mathrm{m}})]}{a-b+bexp[a(z-z_{\mathrm{m}})]}.\end{array}$$(2)

We set the values of the functional parameters following Vangioni et al. (2015): ν = 0.178 M_⊙ ⋅ Mpc⁻³ ⋅ yr⁻¹, the amplitude of the peak of the SFR; z_m = 2, the redshift of the peak of SFR and a = 2.37, b = 1.80, where a is connected to the star formation rate density slope at z < z_m and b − a to the slope for z > z_m. We assume that the fraction of stars in binary systems is 50% (Sana et al. 2012).

While our simulations were performed in metallicity bins, the SFR is given as a function of redshift. In order to calculate the mean metallicity at any given redshift, we used the expression proposed by Belczynski et al. (2016):

$$\begin{array}{c}\hfill \overline{Z}(z)=\frac{y(1-R)}{{\rho }_b}{\displaystyle {\int }_z^{z_{\max }}}\frac{{10}^{0.5}\psi (z^{\prime })}{H_0(1+z^{\prime })\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z^{\prime })}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\mathrm{d}{z}^{\prime },\end{array}$$(3)

where R = 0.27 is the fraction of stellar mass ejected back into the interstellar medium; y = 0.019 is the mass ratio between the new metals created and total stellar mass; and ρ_b = 2.77 ⋅ 10¹¹Ω_bh₀ M_⊙ ⋅ Mpc⁻³ is the baryon density.

At a given redshift, the metallicity content may have a very high dispersion due to differences between galaxies and non-homogeneous mixing within galaxies (we note that we do not use the mass-metallicity relation, as we do not model galaxies by mass). This dispersion influences the environments of stellar formation. For simplicity, we followed Santoliquido et al. (2021) and adopted a log-normal distribution of metallicities around the average metallicity at a given redshift:

$$\begin{array}{c}\hfill P(Z|z)=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp(-\frac{{(log(Z/Z_{\odot })-log(\overline{Z}(z)/Z_{\odot }))}^2}{2{\sigma }^2}),\end{array}$$(4)

with σ = 0.2.

Finally, we could compute the BNS merger rate:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{R}}_{\mathrm{merg}}(t)={\displaystyle \underset{Z_{\min }}{\overset{Z_{\max }}{\int }}}{\displaystyle \underset{t_{\mathrm{delay},\min }}{\overset{t_{\mathrm{delay},\max }}{\int }}}[& \alpha (Z)\psi (t-\mathrm{\Delta }{t}_{\mathrm{delay}})P(\mathrm{\Delta }{t}_{\mathrm{delay}}|Z)\hfill \\ \hfill & P(Z|t-\mathrm{\Delta }{t}_{\mathrm{delay}})\mathrm{d}\mathrm{\Delta }{t}_{\mathrm{delay}}\mathrm{d}Z],\hfill \end{array}\end{array}$$(5)

where the time delay distribution P(Δt_delay|Z) and the fraction of successful BNSs is taken directly from our COSMIC simulations.

3. Results

We show our results below, starting with the three evolutionary tracks we identified using representative systems as examples in Sect. 3.1. We then show in more detail the properties of the progenitor (Sect. 3.2) and BNS populations (Sect. 3.3) for each track, at constant metallicity and across metallicity bins. Finally we use the simulated binaries to obtain a cosmological population and show the resulting merger rate density (Sect. 3.4).

3.1. Evolutionary tracks

3.1.1. Unequal mass ratios: Standard channel

The first evolutionary track we discuss corresponds to Channel I in Vigna-Gómez et al. (2018), Neijssel et al. (2019), Iorio et al. (2023); Channel B in Kruckow et al. (2018, see appendix C) and the dominant channel in Chruslinska et al. (2018), and was also already discussed in the literature as the dominant evolutionary track to produce NSBH systems (e.g. Broekgaarden et al. 2021; Zevin et al. 2020) and BBHs (e.g. Belczynski et al. 2018).

We describe here a representative example of this evolutionary track as shown in Fig. 2. Initially, the binary has an unequal mass ratio (in this example, q ≃ 0.9 is quite on the high end of this distribution) such that the evolution timescale of the primary is shorter than that of the secondary.

Fig. 2. Evolution of a binary system of the standard evolutionary track forming a BNS at Z = 0.00044. The first 24 Myr of main sequence evolution are hidden. The top colour bars show the evolution of stellar types of the primary (top) and the secondary (bottom) during the binary evolution. Then, from top to bottom, we show the evolution of (1) the stellar radii (solid lines) and Roche radii (dashed lines) and (2) the masses of the two stars; (3) the mass growth rates and (4) the semi-major axis (solid line) and the eccentricity (dashed line) of the system. Quantities are shown in green for the primary and in orange for the secondary. When either of the stellar radii exceeds the Roche radius (panel (1)), a simultaneous accretion event occurs on the companion, in opposite colours (panel (3)). Accretion is also possible onto compact objects from stellar winds from the companion. We note that in panel (3), we show the stellar mass growth rate, which accounts for both mass accretion and mass loss simultaneously. When Ṁ is not visible, this means that the star loses mass overall. When Ṁ > 0 and R < RL, we observed accretion of material onto a compact object (NS or WD) from the stellar winds of the companion.

When the primary leaves the MS (after ∼33 Myr), it grows in radius and loses a significant amount of its mass (∼6.5 M_⊙) during a phase of stable mass transfer. This mass is transferred to the secondary. This phase lasts for ≲1 Myr. The primary is left without envelope as a NHeMS, until it again evolves into a NHeHG after ∼5 Myr. At that stage, a brief episode of mass transfer is triggered and the primary explodes in a SN. Given the progenitor mass at the moment of the explosion, the remnant formed is a NS. In this example, after the SN the semi-major axis of the orbit increases (to ∼2500 R_⊙), as well as the orbital eccentricity (e ∼ 0.7). The secondary in turn evolves after a short time (2 Myr) given the mass gained during the first episode of mass transfer. When helium fusion in the core ignites, the secondary enters the giant phase and CE starts ∼2 Myr later, due to the prior eccentricity of the orbit. At the end of the CE phase, the secondary is stripped of its envelope and the semi-major axis is reduced to ∼2.5 R_⊙, while friction during the CE circularises the orbit. Because of the short orbital separation, a continuous case BB mass transfer from the stripped star onto the NS occurs (Tauris et al. 2015). Shortly after, the secondary explodes in an ultra-stripped SN and forms the second NS of the system. Because of the reduced NS kick, the system is given a new initial semi-major axis (∼4 R_⊙) and eccentricity (e ≃ 0.4) but has increased chances of survival. This system eventually merges due to the emission of GWs after ∼5 Gyr. We observe that for systems that follow this track, the stellar evolution time is negligible compared to the typical time for merger due to emission of GWs (see Sect. 3.3).

To summarise, this evolutionary track features a phase of CE between a giant and a NS, and a later phase of mass transfer between the stripped star and the NS (case BB mass transfer). As is discussed in Sect. 4, this track does not produce as many systems as the other evolutionary tracks with our choice of model parameters, in particular at high metallicities.

3.1.2. Equal-mass ratios: Co-evolution

At all metallicities, many BNS systems are produced from the evolution of binary systems whose ZAMS star masses were almost equal (q ≃ 1). In the example shown in Fig. 3, M₁ = 13.44 M_⊙ and M₂ = 13.21 M_⊙, that is q = 0.98. The semi-major axis is initially of ∼2900 R_⊙ and the eccentricity e ≃ 0.69. Both stars evolve on similar timescales and both are in their CHeB phase after 16.44 Myr. Due to the expansion of both stellar envelopes and to the slightly higher mass and thus radius of the primary, a phase of stable mass transfer is initiated, leading to the circularisation of the orbit and therefore reduction of the semi-major axis. Both stars enter in contact and a joint phase of CE occurs, where the two helium cores orbit within the CE. After the envelope ejection, the helium cores have a separation of ∼8.7 R_⊙, and most of the mass has been lost by the binary. The primary is slightly more massive than the secondary, and therefore it evolves first, loses some of its mass to the secondary by case BB mass transfer, and explodes in an ultra-stripped SN. Then the secondary evolves, and stable mass transfer onto the NS follows (case BB mass transfer) until the secondary explodes in another ultra-stripped SN.

Fig. 3. Evolution of a binary system on the equal-mass evolutionary track forming a BNS at Z = 0.0065. The first 14 Myr of main sequence evolution are hidden. See Fig. 2 for panel description.

Several variations of the example described above and shown in Fig. 3 can be found within this evolutionary track, and correspond to nine evolutionary sequences (see Sect. 2.3), but they all share the same distinctive feature of this evolutionary track: a quasi-equal initial mass ratio, leading to a synchronised evolution of both stars and a phase of CE with the two stars on their giant phase. Variations depend on the exact stellar types at the moment of the CE (CHeB or EAGB) and the stellar types of the remaining helium cores at the end of CE. In the cases where the primary is on the EAGB at the onset of CE, the remaining helium cores after CE are already on the HG. A second source of variation is the nature of the stellar type of the helium stars at the moments of the SN. Final variations are found depending on the metallicities, but again, share the same distinct features of this evolutionary track. We note that binary systems in this track typically evolve on a faster timescale than the others, due to the joint evolution of the primary and the secondary. The time between binary formation at ZAMS and BNS formation ranges between ∼10 and 100 Myr.

This evolutionary track is also discussed in the literature, for example in Vigna-Gómez et al. (2018) where it is labelled Channel II; Neijssel et al. (2019) where it is labelled Channel III and Broekgaarden et al. (2021), Iorio et al. (2023) where it is labelled Channel IV. This evolutionary track was also discussed in various contexts by e.g. Brown (1995), Bethe & Brown (1998), Dewi et al. (2006), Hwang et al. (2015). We note that the existence of this evolutionary track is more subject to our choice of model parameters than the others. Indeed, because we impose q_crit = +∞ for mass transfers from stripped stars, this evolutionary track only features one phase of CE. With other choices of parameters, there would be additional variations with one or two additional phases of CE with the helium cores, increasing the number of systems with mergers before the BNS formation and therefore reducing the efficiency of this evolutionary track.

3.1.3. Accretion-induced collapse

The final evolutionary track found with this model features the AIC of a WD into a NS. This evolutionary track has been discussed in Wang & Liu (2020), where the authors show in particular the limited parameter space at ZAMS (P_orb, M₁, M₂) within which this process can occur. Chruslinska et al. (2018) show that this evolutionary track could contribute to the BNS merger rate, owing to lower natal kicks, while Dominik et al. (2012) estimate its contribution to 8% of the total BNS formed at solar metallicity. AIC of WDs is also discussed in a more general context in e.g. Meng et al. (2009), Wang & Han (2010). In our case, this evolutionary track contains two evolutionary sequences as defined in Sect. 2.3. One of them is dominant at the lowest redshifts and differs from the other one only by one element: the primary does not enter the FGB. The physical properties of both evolutionary sequences, however, are similar. We show an example of a representative system in Fig. 4. In this evolutionary track, the initial mass ratio can reach much lower values than in the other tracks, 0.5 ≲ q < 1. In this example, q ≃ 0.48. The initial mass of the secondary is lower than the typical mass of NS progenitors (here M₂ = 4.06 M_⊙), due to the first phase of mass transfer.

Fig. 4. Evolution of a binary system on the AIC evolutionary track forming a BNS at Z = 0.00044. The first 20 Myr of main sequence evolution are hidden. See Fig. 2 for panel description.

The primary evolves first from the MS after 34 Myr. When it reaches the HG, a phase of stable mass transfer starts, reducing the semi-major axis and circularising the orbit. During this event of mass transfer, the primary evolves and loses more mass to the companion. Because the mass transfer is not entirely conservative (some mass is lost by the binary), the semi-major axis increases again, until the primary has lost its entire hydrogen envelope and becomes a stripped helium star. All these events occur on time scales ≲0.1 Myr. At that point, the secondary has a mass that is higher than the initial primary mass. The primary helium star evolves (here in ∼6 Myr), and a new phase of mass loss can occur. When the primary enters the NHeGB, the remaining core mass is therefore much lower than in the other scenarios (here 1.36 M_⊙). At the moment of the SN, the primary collapses into an oxygen-neon WD. As all the mass of the helium core forms the WD, no natal kick is imparted to the WD. The orbit thus remains circular. In turn, the (now massive) secondary leaves the MS and evolves, roughly 20 Myr later. When it reaches the HG, a phase of CE is triggered, with the WD. After the envelope ejection, the orbital separation is reduced to an extremely low value (here 1.3 R_⊙). Throughout the evolution of the helium star, stable mass transfer occurs and gradually increases the mass of the WD. When it reaches the Chandrasekhar mass by accretion, the WD collapses and forms a NS of 1.24 M_⊙ (this value is unique for all systems by construction in COSMIC but would nevertheless be expected to have a low dispersion). Just after, the secondary in turn collapses into a NS thanks to the mass gained during the first phase of mass transfer. The BNS system formed that way has an extremely small initial semi-major axis and merges in typically a few tens of millions of years.

In the context of BBH and NSBH formation history, several authors discuss an additional evolutionary track where no phase of CE is triggered, and all mass transfers are dynamically stable (Neijssel et al. 2019; Broekgaarden et al. 2021; Iorio et al. 2023). In this study, we put our focus on BNS systems and do not find that such an evolutionary track contributes significantly to the population of merging BNSs. This is expected (see e.g. Iorio et al. 2023), since in the case of BBH and NSBH progenitors, their higher masses allow for longer phases of mass transfer, which can shrink the orbit significantly. Additionally, because of the higher remnant masses, systems with higher initial semi-major axes have a lower probability to be disrupted by either of the SNe, and this allows detached binaries, which evolve without a phase of CE, to still produce compact object remnants that remain gravitationally bound.

3.2. Properties of stellar progenitors

We now focus on the properties of the stellar progenitors at fixed metallicity. As could be expected from the discussion in Sect. 3.1, they differ across the three evolutionary tracks. In Fig. 5, we show the distributions of initial binary properties at ZAMS for the metallicities Z = 9.5 × 10⁻⁵, Z = 9.5 × 10⁻⁴ and Z = 1.4 × 10⁻², for each evolutionary track.

Fig. 5. Distribution of the properties of the population of BNS progenitors at their formation (ZAMS) and at three metallicities: Z = 9.5 × 10−5 (left), Z = 9.5 × 10−4 (centre), and Z = 1.4 × 10−2 (right). From top to bottom, the three rows show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes. The standard evolutionary track is shown in blue, the equal-mass track in green, and the AIC track in orange. In black, we show the total population.

In this figure, we filter the simulated populations at each metallicity such that Δt_delay < t_Hubble in order to simplify our comparative studies and to focus only on systems that may merge within the Hubble time. Because binaries with higher metallicities are more likely to form at lower redshift (see Eq. (3)), the tails of these populations with the longest merger times do not contribute to the population of BNSs that merged through cosmic time until today, even if their merger time is lower than the Hubble time. The merger rates (Sect. 3.4) are computed using the exact shapes of these distributions. By construction (see Eq. (5)), they do not include the systems such that t(z_ZAMS)+Δt_delay > t_Hubble, where t(z_ZAMS) is the time at which the progenitor binary is formed.

While the distribution of the various binary properties at ZAMS vary across metallicities, several features are common. We observe in the first and second columns of Fig. 5 that the AIC evolutionary track takes place for quite narrow mass ranges of the primary due to the restrictive conditions for AIC to occur, and slightly wider ranges for the secondary. As expected, these masses are on the lower end, around 5 M_⊙, since these stars are not massive enough to collapse directly into a NS. The equal-mass track occurs primarily for higher-mass progenitors, with the mass distribution extending beyond 20 M_⊙, while the mass distributions for the standard track falls between the equal-mass track and AIC. We also notice that for all tracks the mass distribution shifts to higher values at higher metallicity. This effect is to the first order due to enhanced stellar winds at those metallicities: to form NS progenitors of similar masses right before core-collapse, a higher-metallicity star must have a larger mass at ZAMS than a lower-metallicity star. In principle, overshooting of the convective envelope into the core may also play a role here, but this is modelled in COSMIC as in Pols et al. (1998), without mass or metallicity dependence.

The AIC track disappears almost completely at solar metallicity, as can be seen on the right column of Fig. 5. In fact, as we show in Sect. 3.4 it becomes sub-dominant to the other tracks at z ≲ 3.5. On the other hand, the equal-mass evolutionary track dominates at higher metallicities (see the right column of Fig. 5). To further explore this effect for the AIC track, we show the evolution of the progenitor properties at ZAMS with metallicity in Fig. 6. In fact, at a given metallicity, if M₁ is too high, the primary directly collapses into a NS and if M₁ is too low, it collapses into a carbon-oxygen WD that cannot undergo AIC. The range of masses of oxygen-neon WD progenitors (necessary for AIC events) is therefore quite narrow. When metallicity increases, the mass lost by the primary before its core collapses also increases, thus requiring higher initial ZAMS masses M₁ to produce oxygen-neon WDs. Indeed, as can be seen in Fig. 6, the primary mass increases with metallicity until it stalls around 8 M_⊙. The reason this trend stops is due to the transition between the two evolutionary sequences that compose the AIC evolutionary track (see Sect. 3.1.3). At these higher metallicities, it is this time the secondary whose mass M₂ gradually increases with metallicity. The reason is again related to the mass losses from stellar winds: the secondary initial mass must be higher at higher metallicities. The number of systems produced in the AIC evolutionary track gradually decreases when metallicity increases. This is partially explained by the need for increasing progenitor masses at higher metallicities, which therefore lie further on the decreasing IMF (Kroupa 2001). Combined with the very restrictive mass range of the progenitors of oxygen-neon WDs, and other potential effects linked to binary evolution, this progressively reduces the contribution of the AIC evolutionary track to the total as metallicity increases.

Fig. 6. Distributions of stellar progenitor properties at ZAMS for the AIC evolutionary track and for various metallicities. From top to bottom, the three panels show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes.

The distribution of progenitor properties at ZAMS for the standard and the equal-mass evolutionary tracks are shown in Appendix B. The standard evolutionary track has two components, with similar ZAMS semi-major axes. At lower metallicities, the initial masses M₁ and M₂ are lower and reach up to M₁ ≲ 9 M_⊙; M₂ ≲ 8 M_⊙ at higher metallicities. Above Z ≃ 0.1 Z_⊙, this track almost vanishes. The distribution of progenitor masses for the equal-mass track does not exhibit any particular trend.

Another interesting feature concerning the BNS progenitors is the range of initial semi-major axes. In the first and second columns of Fig. 5, we observe that the AIC track occurs for the lowest progenitor semi-major axes, while the equal-mass track takes places for larger semi-major axes of ∼1000 R_⊙. These findings are in agreement with the progenitor ranges discussed for the AIC track by Wang & Liu (2020), namely M₁ ∼ 6 − 10 M_⊙; q ∼ 0.2 − 0.4 and P ∼ 400 − 1000 days. We also see the impact of metallicity as discussed previously.

Finally, we also note that the initial eccentricity does not affect the tracks followed by the system. The distribution of eccentricities for the progenitors that successfully produce BNSs that merge within the Hubble time is the same as the distribution assumed for the entire population. A possible application of these results would be to find general criteria that would determine which evolutionary track a given binary is more likely to follow (see also the discussion in Broekgaarden et al. 2019). This study is, however, beyond the scope of the present work and we leave it for future investigations.

3.3. Properties of binary neutron star populations

We now discuss the properties of BNS systems at their formation (i.e. right after the secondary SN) for the three evolutionary tracks. We focus on three parameters: Stellar evolution timescale Δt_evol (time from the formation of the stellar binary (ZAMS) to the formation of the BNS), GW orbital decay timescale Δt_GW (time from BNS formation to merger) and the total delay time Δt_delay = Δt_evol + Δt_GW. These parameters are shown in Fig. 7.

Fig. 7. Distribution of the properties of the population of BNS at their formation (i.e. right after the secondary SN) and at three metallicities: Z = 9.5 × 10−5 (left), Z = 9.5 × 10−4 (centre), and Z = 1.4 × 10−2 (right). From top to bottom, the three rows show the distributions of (1) stellar evolution times, (2) merger times, and (3) delay times. As in Fig. 5, the standard evolutionary track is shown in blue, the equal-mass track in green, and the AIC track in orange. In black, we show the total population. The count sum of the three evolutionary tracks accounts for ∼95% of systems (see Sect. 2.3).

From Fig. 7 it is immediately clear that, for individual metallicities, the distribution of delay times does not follow the power-law scaling of 1/t as is often assumed in the literature (e.g. Chruslinska et al. 2018; Neijssel et al. 2019). In fact, depending on the metallicity and on which evolutionary track is dominant, the distribution of delay times can be flat (see the right column of Fig. 7) or extremely peaked (see the left column of Fig. 7, where the counts are shown in logarithmic scale). The equal-mass evolutionary track even shows multiple peaks, which are due to the nine different evolutionary sequences that constitute this track (see Sect. 2.3) and appear for specific values of M₁ at ZAMS. Moreover, a significant fraction of binaries has a delay time dominated by the stellar evolution timescale, with extremely short merger times. Notably, this is the case for the AIC evolutionary track, for which the delay times due to GW emission are extremely short, of the order of 10 Myr, and the total delay time is therefore dominated by the stellar evolution timescale, of the order of 100 Myr. This effect occurs because the distribution of BNS semi-major axes at formation for the AIC track is very peaked around low values, below R_⊙. On the contrary, systems formed via the equal-mass evolutionary track are born with much larger semi-major axes and their delay time is therefore dominated by the orbital decay due to emission of GWs.

Another observation is that the different evolutionary tracks are well distinguished by the progenitor binary evolution time Δt_evol. Progenitors with equal-mass ratios evolve more rapidly due to their initially higher masses (see the green distributions in Fig. 5, first two rows and Fig. 7, first row), BNSs formed through the standard evolutionary track have slightly longer stellar evolution times, while the longest evolution times are required for the AIC evolutionary track, mostly due to the smaller total mass in the system.

The distributions of delay times Δt_delay all have a minimum at ∼1 Myr corresponding to systems with the shortest evolution times and rapid mergers (< 10⁻¹ Myr). Conversely, the maximum delay time is fixed to be t_Hubble, though all these distributions extend to longer delays. The peaks of the distributions for each evolutionary track are primarily affected by the intensities and orientations of the NS natal kicks, but also by events of stellar evolution, such as CE efficiency, case BB mass transfer rates, or the radius expansion of stripped stars. A careful study of the impact of these parameters on the delay time distributions properties would require a systematic variation of those parameters and is outside of the scope of this paper.

To complement the information in Fig. 7, we show the joint distributions of the stellar evolution timescale Δt_evol and the GW orbital decay timescale Δt_GW at a given metallicity Z = 9.5 × 10⁻⁴ in Fig. 8. This joint representation allows one to better visualise the correlation between these two quantities for the three evolutionary tracks. The standard evolutionary track has long evolution times, but not the longest, and a distribution of merger times that ranges between < 1 Myr and t_Hubble. The equal-mass evolutionary track has short evolution times and long merger times, while the AIC evolutionary track has the longest evolution times and short merger times.

Fig. 8. Correlation between the formation and merger times for the three evolutionary tracks at Z = 9.5 × 10−4 and for systems with Δtdelay < tHubble. The three tracks are found on distinct regions of the plane, as discussed in Sect. 3.3. The shaded contours represent the smallest regions containing 10%, 50%, and 90% of the systems.

These results could have profound effects on the properties of BNS populations across cosmic times. Since dominant tracks vary as a function of metallicity, we may expect the BNS properties (in particular their delay time) to also vary with metallicity and hence redshift. In other words, the merging BNS population observed in the local Universe may not be representative of the higher redshift population that merged at an earlier epoch. We further discuss this idea in Sect. 4.

The evolution of BNS properties as a function of progenitor metallicity is presented in Fig. 9. A general trend that we observe is the decrease in the efficiency of BNS formation as metallicity increases. This is mostly due to the increased stellar wind intensities that push progenitors to higher ZAMS masses, further in the tail of the initial mass function distribution. Overall, about one system in 10⁴ − 10⁶ produces a BNS that merges within the age of the Universe.

Fig. 9. Distribution of the properties of the population of BNS at their formation (i.e. right after the secondary SN) for the standard evolutionary track (left), the equal-mass track (centre), and the AIC track (right). Each colour represents the population at one of the sampled metallicity values. From top to bottom, the three rows show the distributions of (1) stellar evolution times, (2) merger times, and (3) delay times.

On the left column in Fig. 9, we observe again that the standard evolutionary track is far less efficient at high metallicities. We confirmed that the standard track produces systems with Δt_delay > t_Hubble, all the more at higher metallicities. These systems are removed in this study. We also observed that the distribution of the initial semi-major axes of the BNSs does not evolve with metallicity; however, the stellar evolution timescale decreases with metallicity. This is expected, since the initial stellar masses of BNS progenitors are higher at higher metallicity, which leads to quicker stellar evolution timescales. We also stress that even in the standard evolutionary track, the distribution of time delays does not follow the 1/t scaling. Indeed, the distribution of initial BNS semi-major axes is not log-normal, and moreover, the stellar evolution timescale is of the same order or longer than the GW timescale, for a significant fraction of the systems.

We now examine the equal-mass track, shown on the central column of Fig. 9. In this case the distributions of timescales are much more complex, in particular the distribution of stellar evolution times. As we mentioned in Sect. 3.1.2, this track is composed of nine evolutionary sequences, which in part leads to the observed complexity, because their relative contributions to the total varies with metallicity. This evolutionary track produces BNS systems efficiently at all metallicities and therefore dominates among the (local) BNS events observed with GW detectors, since the standard and AIC tracks are inefficient at higher metallicities.

The AIC evolutionary track is shown on the right column of Fig. 9. One major difference between this track and the other two is the distribution of merger times. In the AIC track, the semi-major axis after the SN of the secondary is always extremely reduced after the phase of CE with the WD. After the final phases of mass transfer from the helium companion onto the WD that lead to its AIC, and the ultra-stripped SN of the secondary, the semi-major axis still remains extremely small in most cases, leading to a merger in typically ≲10 Myr. This property does not evolve with redshift. We note that the radius evolution of stripped stars is uncertain at low metallicities, and we thus expect these results to be dependent on our choice of model parameters. Only the overall number of BNSs formed by this evolutionary track decreases at higher metallicities. The standard evolutionary track has much longer merger times at low metallicities, while for the equal-mass evolutionary track, the merger time tends to increase with metallicity.

3.4. Merger rate density

We present in Fig. 10 the merger rate obtained after combining all the metallicity bins following the procedure described in Sect. 2.4, along with the constraint on the local merger rate from GW observations: ℛ_merg(z = 0) = 10 − 1700 Gpc⁻³ ⋅ yr⁻¹ (Abbott et al. 2023b). Our model is compatible with the observed local merger rate constraint, although we note that it still covers a broad range. In that regard, new detections will help constrain population models.

We now focus on the contributions from each of the three evolutionary tracks. First, we observe that the standard evolutionary track contributes to a very small fraction of the merger rate at all redshifts. Moreover, the evolutionary track dominating the population of merging binaries changes with redshift: at higher redshifts, most mergers are produced by BNS systems formed via the AIC evolutionary track, while at lower redshifts, the population we see merging in the local Universe is mostly formed by equal-mass ZAMS progenitors. This result suggests that the BNS population observable with GWs may not be representative of the entire merging BNS population across cosmic times.

The two peaks in the merger rate contribution from the AIC evolutionary track originate from the two evolutionary sequences that contribute to this track that we discussed in Sect. 3.1.3: one at low metallicities and the other at higher metallicities. In the transition between the two peaks, the evolutionary sequence where the primary goes through the FGB (higher redshifts) is not as efficient, and the second sequence where the primary does not enter the FGB phase is not yet efficient.

Because most merging systems at high redshift are formed by the AIC evolutionary track, they have very short merger times Δt_GW. The merger rate density shown in Fig. 10 therefore closely follows the evolution of the SFR density.

Fig. 10. Evolution of the merger rate with redshift. We show the total rate (black), standard evolutionary track (blue, see Sect. 3.1.1), equal-mass progenitors (green, Sect. 3.1.2), and AIC (orange, see Sect. 3.1.3). The relative contribution of each track to the merging population varies with redshift, with the AIC track prominent at high redshifts and the equal-mass track dominant at low redshifts. The grey box represents the constraint on the merger rate in the local Universe after the O3 observing run of LIGO/Virgo (Abbott et al. 2023b). It is stretched up to z = 0.3 for visualisation purposes.

Finally, it is interesting to note that the total merger rate evolution is quite smooth and hides the diversity of progenitor evolutionary tracks that make up the BNS population, and despite the broad diversity of the distribution shapes shown in Fig. 7.

4. Discussion

We have explored the formation scenarios of merging BNSs using the population synthesis code COSMIC for one chosen set of physical parameters (Sect. 2). We have discussed the properties of the stellar progenitors at different metallicities as well as the resulting BNS populations and studied the predicted merger rates. Our results are summarised as follows:

• We identified three dominant evolutionary tracks: the standard track, also commonly studied in other works, that involves a phase of CE with the first-born NS followed by an episode of case BB mass transfer; the equal-mass track, where both progenitor stars evolve on synchronised timescales and which includes a phase of CE when both stars are in their giant phase; and a track that involves the formation of an O-Ne WD followed by its AIC due to mass accretion from the companion.

• With our choice of model, we found that the equal-mass evolutionary track is the major contributor to the population of merging BNSs at low redshift (higher metallicities), while the AIC track is dominant at high redshift (lower metallicities). The standard track, in contrast, is always sub-dominant, as it mostly forms systems with Δt_delay > t_Hubble.

• Similar to previous studies, we find that BNSs that merge within the Hubble time undergo at least one episode of CE. We also find a merger rate at z = 0, which is in agreement with the constraints from current LIGO/Virgo observations, although more detections are needed to obtain better constraints. The merger rate peaks at redshift z ≈ 1.7.

• The distribution of time delays between the formation of stellar progenitors and the BNS merger does not follow a simple power law at a fixed progenitor metallicity. Instead, the distribution is typically complex, its shape evolves with metallicity and depends on the relative contribution of the different evolutionary tracks.

• In some cases, particularly at low metallicities, the stellar evolution time (the time to form a BNS) can be longer than the orbital decay time of the BNS due to emission of GWs.

There are two aspects of our results that deviate from previous similar studies of BNS populations, namely Kruckow et al. (2018), Belczynski et al. (2018), Vigna-Gómez et al. (2018), and Iorio et al. (2023), which used the COMBINE, StarTrack, COMPAS, and SEVN population synthesis codes, respectively. These studies find that the standard track is dominant at all redshifts, and they do not exhibit the AIC track. One reason for these discrepancies may be the difference of population synthesis codes used in these studies as well as different choices of physical parameters. COMPAS, used by Vigna-Gómez et al. (2018), is the closest population synthesis code to COSMIC since it is also based on BSE. However, Vigna-Gómez et al. (2018) focused on Galactic systems without selecting merging binaries, whereas in this work we only discuss systems merging within the age of the Universe. In view of our results, it is plausible that these two distinct populations could have different formation channels, and we plan to explore this topic further in future work. In particular, the AIC track that only appears at high redshift produces very short merger-time systems and may therefore be exclusively related to merging systems. Belczynski et al. (2018) focused on merging systems and assumed that all the progenitor stars formed at a single metallicity of Z = 0.01 in their isolated binary formation scenario. The dominant evolutionary track in their work corresponds to our standard track, with two phases of CE owing to an unstable mass transfer from the stripped helium companion with the first-born NS. In this work, we find a dominant contribution from the equal-mass track at such a slightly sub-solar metallicity. This difference could be due to a different choice of binary evolution parameters and will be explored in a follow-up study. Chruslinska et al. (2018) show that the AIC evolutionary track may contribute to the population of merging BNSs, though it is sub-dominant at the three metallicities they sample. Similarly, Dominik et al. (2012) show a contribution of the AIC track of 8% at solar metallicity, which disappears at 0.1 Z_⊙. These results are also at odds with our findings. As discussed in Ruiter et al. (2019), the properties of BNS systems formed in the AIC track may be particularly sensitive to CE physics, but the contribution of this track may be significant. Finally, we stress that COSMIC relies on pre-computed single stellar evolution tracks. Several recent population synthesis codes feature their own evolutionary tracks, most notably POSYDON (Fragos et al. 2023) and SEVN (Iorio et al. 2023), which are based on more recent results using MESA (Jermyn et al. 2023 and references therein) and PARSEC (Bressan et al. 2012). They rely on interpolations between the tracks rather than on fitting formulae and lead to substantial differences in the populations of formed compact objects. In light of these results, a dedicated comparison between the various population synthesis codes would be needed in order to fully understand the sources of the discrepancies we mentioned.

While we did not perform a full parameter study, we looked into different values for the critical mass ratios q_crit, using qcflag (see Sect. 2). We observed that this flag has a strong effect on binary evolution. When setting qcflag = 1 or 2, the mass transfers from stripped helium stars are often dynamically unstable (no case BB mass transfer) and lead to a new phase of CE, while in the model presented here (qcflag = 5), they are always dynamically stable. Our choice of model therefore naturally produces evolutionary tracks with only one CE phase, while the others have two, even sometimes three, CE events.

We also tested another aspect of CE physics. In our model, we assumed that stars automatically merge during CE if the companion has no core-envelope boundary (cemergeflag = 1). We tested the effect of this assumption by allowing such binaries to survive CE (cemergeflag = 0), and we found that this typically adds a few additional evolutionary sequences. In particular, in this alternative scenario, there is an additional AIC evolutionary sequence that dominates at high metallicities. A CE then appears when the primary has already formed a WD and the secondary is on the HG. These preliminary tests emphasise the importance of a full-scale parameter study, which we leave for future work.

In this work, we do not address the question of NS masses. Indeed, there are only two observed merging BNSs as of writing. The observed Galactic population of NSs gives more insight into their actual mass distribution, and if coupled with another treatment of NS mass in population synthesis codes (as in, e.g., COMBINE; Kruckow et al. 2018), it could provide additional parameter constraints. This interesting aspect should be investigated more in the future.

The results presented in this work, if confirmed, could be important for the study of the r-process abundances in the interstellar medium. Indeed, the multi-messenger observations of GW 170817 have proved that r-process elements are synthesised in the ejecta of merging BNSs and power kilonova emission (Abbott et al. 2017b). It seems natural to assume that all the r-process elements in the Galaxy were forged in BNS mergers since the overall merger rates are consistent with the estimated total mass in r-process elements, though the predictions for the exact yields per event are still uncertain, owing to the effect of the equation of state or the ejecta geometry, for example. However, detailed semi-analytic and hydrodynamical models of the Milky Way show that in order for BNS mergers to be the dominant source of r-process elements in metal-poor stars, the delay time between the formation of progenitor stars and BNS merger Δt_delay should be very short, on the order of a few, or at most, tens of millions of years (e.g. Beniamini et al. 2018; Beniamini & Piran 2019; van de Voort et al. 2020; Kobayashi et al. 2023). Such delays are much shorter than those predicted for the standard evolutionary track. In our study, BNSs formed via the AIC track have very short delay times and are also more commonly formed at lower metallicities. Because the AIC events do not eject final fusion products such as iron, unlike core-collapse SNe, this evolutionary track allows for the local enrichment of the surrounding gas in r-process elements, with a limited enrichment in nuclear fusion products. This could therefore naturally explain the existence of low-metallicity stars enriched with r-process elements.

Statistical studies of the offsets of short gamma-ray bursts (GRBs) in their host galaxies (e.g. Fong et al. 2022) show that a fraction of short GRBs are observed close to the centre of their host galaxies. Such observations are indirect evidence that progenitors of short GRBs (often associated with BNS mergers) have low natal kicks or short merger times. Other studies have pointed out that most short GRB host galaxies are still star forming (Nugent et al. 2022), which again hints towards a population of short delay time BNS mergers. Our results naturally provide such a population, mostly in the AIC track.

Finally, our results could be important for the next generation of GW detectors. Indeed, while current interferometers can only observe BNSs in the local Universe, the planned third-generation detectors, Einstein Telescope and Cosmic Explorer (Maggiore et al. 2020; Evans et al. 2021; Branchesi et al. 2023), will see much farther out, reaching the cosmic noon and beyond. These observations will therefore give access to a much wider range of environments of BNS progenitors and help uncover their origins.

Acknowledgments

CP acknowledges funding support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université. This work made use of the Infinity cluster, hosted by the Institut d’Astrophysique de Paris, on which the simulations were run and post-processed. The authors warmly thank S. Rouberol for running it smoothly. The authors thank K. Breivik for useful discussions and the entire COSMIC team for making the code publicly available.

References

Appendix A: Parameters used in the COSMIC simulations

We list in Table A.1 the parameter values used in our simulation (see Sec. 2 for a more detailed discussion on our choice of parameters). The version of COSMIC used in this work is v. 3.4.0.

Appendix B: Evolution of ZAMS progenitor properties with metallicity

We show here the distributions of the progenitor properties at ZAMS for the standard evolutionary track (Fig. B.1) and the equal-mass evolutionary track (Fig. B.2). These figures are discussed in Sect. 3.2.

Fig. B.1. Distributions of stellar progenitor properties at ZAMS for the standard evolutionary track and for various metallicities. From top to bottom, the three panels show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes.

Fig. B.2. Same as Fig. B.1 but for the equal-mass evolutionary track.

All Tables

All Figures

Fig. 1. Illustration of our classification scheme. Shown is the number of BNS systems for each sequence of stellar types of the primary (rows) and of the secondary (columns) at Z = Z⊙ = 0.014. The combination of primary and secondary evolutionary sequences creates 46 evolutionary sequences, but 95% of the systems are found in 12 of them. The sequences of numbers on the y axis and on the x axis are the series of stellar types (see Table 1) for the primary and the secondary, respectively. Only systems with Δtdelay < tHubble are shown.

In the text

Fig. 2. Evolution of a binary system of the standard evolutionary track forming a BNS at Z = 0.00044. The first 24 Myr of main sequence evolution are hidden. The top colour bars show the evolution of stellar types of the primary (top) and the secondary (bottom) during the binary evolution. Then, from top to bottom, we show the evolution of (1) the stellar radii (solid lines) and Roche radii (dashed lines) and (2) the masses of the two stars; (3) the mass growth rates and (4) the semi-major axis (solid line) and the eccentricity (dashed line) of the system. Quantities are shown in green for the primary and in orange for the secondary. When either of the stellar radii exceeds the Roche radius (panel (1)), a simultaneous accretion event occurs on the companion, in opposite colours (panel (3)). Accretion is also possible onto compact objects from stellar winds from the companion. We note that in panel (3), we show the stellar mass growth rate, which accounts for both mass accretion and mass loss simultaneously. When Ṁ is not visible, this means that the star loses mass overall. When Ṁ > 0 and R < RL, we observed accretion of material onto a compact object (NS or WD) from the stellar winds of the companion.

In the text

	Fig. 3. Evolution of a binary system on the equal-mass evolutionary track forming a BNS at Z = 0.0065. The first 14 Myr of main sequence evolution are hidden. See Fig. 2 for panel description.
In the text

	Fig. 4. Evolution of a binary system on the AIC evolutionary track forming a BNS at Z = 0.00044. The first 20 Myr of main sequence evolution are hidden. See Fig. 2 for panel description.
In the text

Fig. 5. Distribution of the properties of the population of BNS progenitors at their formation (ZAMS) and at three metallicities: Z = 9.5 × 10−5 (left), Z = 9.5 × 10−4 (centre), and Z = 1.4 × 10−2 (right). From top to bottom, the three rows show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes. The standard evolutionary track is shown in blue, the equal-mass track in green, and the AIC track in orange. In black, we show the total population.

In the text

	Fig. 6. Distributions of stellar progenitor properties at ZAMS for the AIC evolutionary track and for various metallicities. From top to bottom, the three panels show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes.
In the text

Fig. 7. Distribution of the properties of the population of BNS at their formation (i.e. right after the secondary SN) and at three metallicities: Z = 9.5 × 10−5 (left), Z = 9.5 × 10−4 (centre), and Z = 1.4 × 10−2 (right). From top to bottom, the three rows show the distributions of (1) stellar evolution times, (2) merger times, and (3) delay times. As in Fig. 5, the standard evolutionary track is shown in blue, the equal-mass track in green, and the AIC track in orange. In black, we show the total population. The count sum of the three evolutionary tracks accounts for ∼95% of systems (see Sect. 2.3).

In the text

	Fig. 8. Correlation between the formation and merger times for the three evolutionary tracks at Z = 9.5 × 10−4 and for systems with Δtdelay < tHubble. The three tracks are found on distinct regions of the plane, as discussed in Sect. 3.3. The shaded contours represent the smallest regions containing 10%, 50%, and 90% of the systems.
In the text

Fig. 9. Distribution of the properties of the population of BNS at their formation (i.e. right after the secondary SN) for the standard evolutionary track (left), the equal-mass track (centre), and the AIC track (right). Each colour represents the population at one of the sampled metallicity values. From top to bottom, the three rows show the distributions of (1) stellar evolution times, (2) merger times, and (3) delay times.

In the text

Fig. 10. Evolution of the merger rate with redshift. We show the total rate (black), standard evolutionary track (blue, see Sect. 3.1.1), equal-mass progenitors (green, Sect. 3.1.2), and AIC (orange, see Sect. 3.1.3). The relative contribution of each track to the merging population varies with redshift, with the AIC track prominent at high redshifts and the equal-mass track dominant at low redshifts. The grey box represents the constraint on the merger rate in the local Universe after the O3 observing run of LIGO/Virgo (Abbott et al. 2023b). It is stretched up to z = 0.3 for visualisation purposes.

In the text

	Fig. B.1. Distributions of stellar progenitor properties at ZAMS for the standard evolutionary track and for various metallicities. From top to bottom, the three panels show the distributions of (1) primary masses, (2) secondary masses, and (3) initial semi-major axes.
In the text

	Fig. B.2. Same as Fig. B.1 but for the equal-mass evolutionary track.
In the text