Double Neutron Star Mergers and Fast Blue Optical Transients: A Circumbinary Disk Connection

INTRODUCTION

Binaries are common across astrophysical scales, ranging from galaxy pairs to stellar and sub-stellar systems. In gas-rich environments, binarity often leads to the formation of circumbinary disks (CBDs), in which material settles into orbit around the binary as a whole rather than around either component individually (e.g., [30, 3, 56, 110, 81, 38]). The exchange of angular momentum and mass between the CBD and the binary through gravitational and accretion forces affects the binary’s evolution (e.g., [5] (Artymowicz & Lubow 1996)), driving either orbital decay or expansion.

The formation of CBDs around stellar binaries and their effect on the subsequent evolution of orbital separation and eccentricity has been extensively studied in recent years, with growing interest in CBDs around post-common-envelope evolution (CEE) binary systems (e.g., [66, 24, 15, 61, 37] (Vetter et al. 2025)). If a binary survives dynamically unstable mass transfer, a fraction of the envelope material may remain bound to the system and fall back towards the inner binary. Since this gas carries substantial angular momentum, it cannot be re-accreted directly onto either component. Instead, it circularizes at radii comparable to a few times the post-CEE orbital separation, forming a rotationally supported CBD

Corresponding author: Aldana Grichener

[email protected] that can persist while it viscously spreads and is slowly accreted (e.g., see Pringle [109]; Lai & Muñoz [76] for reviews).

For binaries that involve a neutron star (NS) or a black hole (BH), the interaction with a CBD formed around the compact object and exposed core of its companion after CEE (see Figure 1) could affect the formation and merger rates of double compact objects (e.g., Vetter et al. [150]; Wei et al. [160]; Unger et al. [145]). If angular momentum is transferred from the binary to the disk, the orbit shrinks, potentially bringing its components into the gravitational-wave (GW)–driven merger regime, or shortening the merger timescale in systems that were already supposed to merge. If the compact object and core of the giant get even closer due to their interaction with the CBD, they may ultimately merge, producing energetic transients such as luminous fast blue optical transients (LFBOTs; e.g., Soker 2022; Metzger [92]; Tuna & Metzger [143]; Cohen & Soker [19]) and Ibn/Icn supernovae (SNe; e.g., Wei et al. [160]); see Grichener [48] for a recent review. Conversely, if the binary system extracts angular momentum from the disk, the orbit widens, driving the two components further apart. This orbital widening can delay or even prevent the merger of double compact objects that would otherwise coalesce within the age of the Universe through GW emission.

The effect of CBDs formed at the end of the CEE phase on binary compact object formation and mergers is particularly important for double neutron stars (DNSs). Most DNSs form through isolated binary evolution following a CEE phase between a giant star and an NS (e.g., Tutukov & Yungelson [144]; Dewi & van den Heuvel [25]), after the exposed core of the giant explodes in a second SN event (e.g., Tauris et al. [131]; Andrews et al. [4]; Tauris et al. [132]; Vigna-Gómez et al. [153]; Fragos et al. [34]; Vigna-Gómez et al. [154]; Laplace et al. [78]; Vartanyan et al. [149]; Grichener [47]). Altering their merger timescales and rates due to CBD interactions could have important implications for future GW observations (The LIGO Scientific Collaboration et al. 2025), observations of Galactic DNSs through pulsar timing (Manchester et al. [86]), and the expected distribution of short gamma-ray bursts (Goldstein et al. [45]; Mandel et al. [89]) and related transients, as well as for galactic chemical evolution (e.g., Kobayashi et al. [70]).

In this work, we incorporate the formation of CBDs right after the CEE phase, and their influence on the orbital evolution of the inner binaries, into the population synthesis code COMPAS (Stevenson et al. [129]; Vigna-Gómez et al. [153]; Team COMPAS: J. Riley et al. [134]; Team COMPAS: I. Mandel et al. [133]). Then, we study the effect of these CBDs on the evolution and mergers of DNS progenitors. In Section 2, we describe our CBD framework and population models. We then explore the implications of CBD interactions on DNS mergers (Section 3) and the occurrence of LFBOTs through ‘premature’ mergers (Section 4). In Section 5, we examine the sensitivity of our findings to different population synthesis assumptions. We find the resulting DNS chirp-mass distributions which may break degeneracies between model assumptions in Section 6. We discuss the implications of the reduction in prompt DNS mergers for r-process enrichment and the caveats of our work in Section 7 and summarize in Section 8.

Methods

CBD Framework in Population Synthesis

We build a formalism to study the effect of CBDs formed around the central binary from leftover envelope that remains bound after the CEE phase on a population level, and integrate it into the rapid population synthesis code COMPAS (Stevenson et al. [129]; Vigna-Gómez et al. [153]; Team COMPAS: J. Riley et al. [134]; Team COMPAS: I. Mandel et al. [133]). To model the evolution of single stars, COMPAS uses the analytic fits from Hurley et al. [59] (see Tout et al. [142]), based on the stellar models of Pols et al. [108]. Binary interactions are modeled semi-analytically following the general approach of Hurley et al. [60] (see also Tout et al. [141]), but with significant changes in the treatment of mass transfer, CEE, winds, SNe, etc., as summarized by Team COMPAS: J. Riley et al. [134]; Team COMPAS: I. Mandel et al. [133]. We incorporate the effect of CBDs by combining an angular momentum based criterion for CBD formation ([145]) and the results of hydrodynamical simulations that study the interaction of the disks with the binaries they surround ([120 b]; as interpolated by Valli et al. [146]). Our framework is illustrated in Figure 2.

Vector (PDF)

To determine whether a system can form a CBD, we evaluate $j_{\mathrm{avail}}$, the post-CEE specific angular momentum retained by the bound envelope material and available for disk formation. It depends on the fraction $f_M$ of the pre-Roche-lobe-overflow (RLOF) giant’s envelope mass $M_{\mathrm{env}}$¹ that remains bound and may form a CBD

of mass

and on the fraction $f_J$ of the total angular momentum of the binary $J_{\mathrm{total}}$ at the onset of RLOF that this mass retains

Under this parametrization, $j_{\mathrm{avail}}$ is given by

where $\beta\equiv f_J/f_M$ is a parameter that quantifies the specific angular momentum retention of the bound envelope material.

We compare $j_{\mathrm{avail}}$ to the Keplerian specific angular momentum required to support a CBD

where $M_{1,\mathrm{post-CEE}}$ and $M_{2,\mathrm{post-CEE}}$ are the masses of the binary components that survived the CEE phase, $R_{\mathrm{CBD,in}}$ is the inner disk radius of the potential CBD, and $\zeta$ is a dimensionless parameter that accounts for the mass distribution of the disk. We take $R_{\mathrm{CBD,in}} = 2.5a_{\mathrm{post,CEE}}$, where $a_{\mathrm{post,CEE}}$ is the post-CEE orbital separation ((e.g., Kashi & Soker 2011; Tuna & Metzger 2023; inspired by Artymowicz & Lubow 1994)), and adopt $\zeta= 4$, which is roughly consistent with the density profile obtained in the CBD simulations of Vetter et al. [150]. If

we assume that the leftover envelope material has enough specific angular momentum to form a centrifugally supported CBD around the post-CEE binary.

For systems where the condition in Equation (5) is satisfied, we incorporate the effect of CBDs on the evolution using SPINDLER (Valli et al. [146]), an interpolation across the hydrodynamical simulations of Siwek et al. [120], and additionally apply corrections that account for super-Eddington accretion. Siwek et al. [120] perform a large suite of two-dimensional scale-free hydrodynamical simulations of binaries immersed in CBDs using the moving-mesh code AREPO ([128, 95, 103]). The disk is treated using an $\alpha_{\mathrm{visc}}$-viscosity prescription ([118]) with $\alpha_{\mathrm{visc}} = 0.1$ and fixed aspect ratio $h = 0.1$. The disk mass is assumed to be small compared to the binary mass, so gravitational instability is not expected (though see discussion in § 7.2.1). The binary components are represented by sink particles and kept on fixed Keplerian orbits, allowing the simulations to measure not only the gravitational and accretion torques exerted by the gas flow, but also the accretion onto each component individually. The secular evolution of the binary is then reconstructed from these quantities. Valli et al. [146] interpolate across the results of these hydrodynamical simulations, and develop a model that predicts the orbital separation, eccentricity and mass ratio of a given binary following its interaction with a CBD.

The parameter $\beta$ affects the interaction only through the CBD formation criterion: large values of $\beta$ correspond to greater retained specific angular momentum and therefore allow CBDs to form in a larger fraction of systems. By contrast $f_{\mathrm{M}}$ sets the mass of the CBD, and therefore the amount of material available to accrete onto and exchange angular momentum with the binary. The inner binary in our simulations can accrete a fraction $f_{\mathrm{acc}} \leq f_{\mathrm{M}}$ from the pre-CEE envelope mass

We use SPINDLER to determine how the total accreted mass $M_{\mathrm{acc}}$ is partitioned between the two binary components, and subsequently correct the accreted amounts according to each stellar type. For the exposed core, accretion is governed mainly by disk viscosity, and its relatively long thermal timescale allows most of the inflowing material to be accreted. If one of the stars in the inner binary is a compact object, however, super-Eddington accretion is expected to drive winds from an accretion disk around the compact object fed by the CBD (also known as ‘minidisk’, see Figure 1), expelling some fraction of the inflowing material. To check whether the SPINDLER supplied mass corresponds to a super-Eddington accretion rate, we assume a constant mass accretion rate over the CBD lifetime, $\tau_{\mathrm{CBD}}$. We take a characteristic disk lifetime of $\tau_{\mathrm{CBD}} = 10^{5}\ \mathrm{yr}$, which corresponds to a timescale in which the outer edge of the disk has time to viscously spread to radii where photoevaporation becomes significant, and thereby drives wind mass loss that dissipates the disk (e.g., [143, 61]). The CBD interaction proceeds in parallel with stellar and binary evolution, including wind mass loss and mass transfer, until the disk is depleted or the core undergoes a SN, which is assumed to eject the remaining CBD. The assumed timescales imply, however, that in most cases the CBD would be depleted well before the stripped core explodes in a stripped-envelope SN, forming the second NS in the DNS binary.

We compare the constant accretion rate inferred from the SPINDLER supplied mass with the compact object’s Eddington accretion rate. In case of super-Eddington accretion, we consider two prescriptions for the resulting compact-object accretion rate. First, we follow the results of Ghodla & Eldridge [42], which allow for a strong super-Eddington accretion of the form

where $\dot{M}_{\mathrm{acc,CO}}$ is the actual accretion rate by the compact object after accounting for mass loss through super-Eddington winds, $\dot{M}_{\mathrm{supplied,CO}}$ is the SPINDLER-inferred accretion rate, and $f_{\mathrm{GE}}$ is the suppression factor, which is equal to $2/3$ in the case of NSs. To explore a more pessimistic scenario, we take the prescription of Tuna & Metzger [143] which results in

where $\dot{M}_{\mathrm{Edd}}$ is the Eddington mass accretion rate. In both cases, since the mass is lost from the ‘minidisk’ surrounding the accreting compact object, we assume that it carries away the accretor’s specific orbital angular momentum, and we remove the corresponding angular momentum from the binary together with the lost mass. Throughout the mass accretion and loss process, we update the post-CEE orbital separation, eccentricity, mass ratio, and component masses using the values calculated with SPINDLER in each timestep.

Population Synthesis Models

After implementing the CBD formalism presented in § 2.1, we use COMPAS to follow the evolution of massive binaries and study the effect of CBD interaction on the evolution and merger of DNS progenitors. Our ‘fiducial no CBD’ model (model 1) follows the COMPAS code defaults (Team COMPAS: I. Mandel et al. [133]), as described in more detail below. We run a grid of populations varying different evolutionary assumptions (models 2-9), and build a ‘fiducial CBD’ model by adding CBD interactions according to the formalism described in § 2.1 to the ‘fiducial no CBD’ case (model 10). Given the uncertainties in the retained specific angular momentum and the CBD mass, we evolve binaries across a range of CBD parameter choices (models 11-47). Then, we incorporate CBD interactions into the variations explored in models 2-9, allowing us to assess CBD effects across the different stellar and binary evolution assumptions (models 48-55). Each population contains $10^{7}$ massive binaries. The full list of simulations is described in Table 1.

Number	Physics	Variation
1	‘Fiducial no CBD’
2	Metallicity distribution	$f_Z(Z)$ = log − uniform
3	SN recipe	Fryer et al. [35] remnant masses, bimodal kick
4	Mass transfer stability	$q_{\rm crit}$ from Claeys et al. [18]
5	Mass transfer stability	‘adiabatic’ $q_{\rm crit}$ from Ge et al. (2020)
6	Mass transfer stability	‘isentropic’ $q_{\rm crit}$ from Ge et al. (2020)
7	CEE	$\alpha= 0.5$
8	CEE	$\alpha= 2$
9	CEE	$e_{\rm post,CEE} = 0$
10	‘Fiducial CBD’	$\beta= 3$, $f_{\rm M} = 0.05$
11–22	CBD disk mass fraction	$f_{\rm M} = \{0.005, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.095, 0.10\}$, $\beta= 3$
23–46	CBD specific angular momentum	$\beta= \{0.5, 1, 1.25, 1.4, 1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.9, 2, 2.15, 2.25, 2.4, 2.5, 2.75, 3, 3.5, 4, 4.5, 5, 5.5, 6\}$, $f_{\rm M} = 0.05$
47	Super-Eddington accretion mode	Tuna & Metzger [143] super-Eddington accretion, $\beta= 3$, $f_{\rm M} = 0.05$
48	Metallicity distribution, CBD	$f_Z(Z)$ = log − uniform, $\beta= 3$, $f_{\rm M} = 0.05$
49	SN recipe, CBD	Fryer et al. [35] remnant masses, bimodal kick, $\beta= 3$, $f_{\rm M} = 0.05$
50	Mass transfer stability, CBD	$q_{\rm crit}$ from Claeys et al. [18], $\beta= 3$, $f_{\rm M} = 0.05$
51	Mass transfer stability, CBD	‘adiabatic’ $q_{\rm crit}$ from Ge et al. (2020), $\beta= 3$, $f_{\rm M} = 0.05$
52	Mass transfer stability, CBD	‘isentropic’ $q_{\rm crit}$ from Ge et al. (2020), $\beta= 3$, $f_{\rm M} = 0.05$
53	CEE, CBD	$\alpha= 0.5$, $\beta= 3$, $f_{\rm M} = 0.05$
54	CEE, CBD	$\alpha= 2$, $\beta= 3$, $f_{\rm M} = 0.05$
55	CEE, CBD	$e_{\rm post,CEE} = 0$, $\beta= 3$, $f_{\rm M} = 0.05$

To focus on massive stars that could explode as SNe forming an NS remnant, we sample the zero-age main sequence (ZAMS) mass of the primary (initially heavier) star from the Kroupa initial mass function (IMF) $dN/dM_{1,\mathrm{ZAMS}} \propto M_{1,\mathrm{ZAMS}}^{-2.3}$ (Kroupa [74]) with masses between $5 \le M_{1,\mathrm{ZAMS}}/M_\odot\le100$. We draw the mass of the secondary (initially lighter star) from a flat distribution in mass ratio with $0.1 \le q_{\mathrm{ZAMS}} \equiv M_{2,\mathrm{ZAMS}}/M_{1,\mathrm{ZAMS}} \le1$ (Sana et al. [115]) and the initial orbital separation from a flat-in-the-log distribution with $0.1 \le a_{\mathrm{ZAMS}}/\mathrm{au} \le1000$ (Opik 1924). For simplicity, we assume circular orbits at the beginning of the evolution (see de Mink & Belczynski [23] for a study showing the negligible impact of initial eccentricities on the merger rates of double compact objects).

All of the stars in our fiducial model have solar metallicity of $Z_\odot= 0.0142$ (Asplund et al. [7], see discussion in § 7.2.2). We relax this assumption by running models with a log-uniform metallicity distribution with a minimum value of $Z_{\min} = 10^{-4}$ and a maximum value of $Z_{\max} = 0.03$ (e.g., Neijssel et al. [97]; models 2, 48). This distribution is dominated by low-$Z$ systems, as it heavily weights subsolar metallicities.

NSs are formed in SN explosions that eject mass and may impart a natal kick to the compact remnant. In our fiducial models, we use the Mandel & Müller [87] prescription in which the natal kick amplitude depends on the pre-SN carbon-oxygen and helium shell mass. We run an additional model using the Fryer et al. [35] delayed SN prescription and adopting a bimodal kick velocity distribution. We take a higher mode of $\sigma_{\rm high} = 217\ \mathrm{km\ s^{-1}}$ for standard core-collapse SNe (one of the distributions from Disberg & Mandel [26]; correction to Hobbs et al. [57]) and a lower mode $\sigma_{\rm low} = 30\ \mathrm{km\ s^{-1}}$ for electron capture SNe and for ultra stripped SNe where $\lesssim0.2M_\odot$ of the envelope mass is bound to the progenitor at the moment of explosion (e.g., Vigna-Gómez et al. [153]; models 3, 49).

The vast majority of binaries that evolve towards DNS systems go through a CEE phase between the giant secondary and the NS that resulted from the SN explosion of the primary. The CEE phase is initiated when mass transfer from the giant star to its NS companion becomes dynamically unstable. We determine the stability of the mass transfer episode by comparing the donor’s adiabatic response to mass loss with that of the Roche-lobe radius $R_{\rm RL}$. If the donor’s response, quantified by $d\ln R_{\rm donor,MT}/d\ln M_{\rm donor,MT}$, where $R_{\rm donor,MT}$ and $M_{\rm donor,MT}$ are the radius and mass of the donor at the onset of mass transfer, is smaller than that of the Roche lobe, $d\ln R_{\rm RL}/d\ln M_{\rm donor,MT}$, the mass transfer episode is unstable and a CEE phase begins. We perform additional simulations varying the criteria of mass transfer stability, using critical mass ratios that separate between stable and unstable RLOF from Claeys et al. [18] (models 4, 50) and Ge et al. 2020 (models 5-6, 51-52).

To determine the evolution of the binary during the CEE phase, we use the energy formalism (e.g., Webbink [158]; Livio & Soker [83]), which compares the binding energy of the envelope prior to the CEE phase with the difference in orbital energies at the onset and at the termination of CEE, $E_{\mathrm{bind}} = \alpha\Delta E_{\mathrm{orb}}$. In our fiducial model, the efficiency parameter $\alpha$, which represents the fraction of the orbital energy budget available to unbind the envelope, is equal to one. We run two additional simulations with $\alpha= 0.5$ (models 7, 53) and $\alpha= 2$ (models 8, 54). The latter represents an additional energy source needed to unbind the envelope, and was favored by several double compact object studies (e.g., Fragos et al. [34]; García et al. [40]; Zevin et al. [162]; Broekgaarden & Berger [12]). We adopt the standard COMPAS prescriptions for CEE outcomes, without accounting for the effect of CBD torques on mass transfer (see discussion in § 7.2.2).

CEE inspiral typically tends to circularize the orbit through strong dissipation and drag, leaving at most a small residual eccentricity in most surviving post-CE binaries (e.g., Ricker & Taam [113]; Sand et al. [116]; González-Bolívar et al. [46]; see discussion in § 7.2.2). Due to the uncertainty in the exact eccentricity and its potential importance to the interaction of the CBD with the binary, in the fiducial case we adopt a simple phenomenological prescription in which the post-CEE eccentricity is set to 10% of its pre-RLOF value. We relax this assumption and run additional models with circular post-CEE binaries (models 9, 55), following standard practice in population synthesis studies.

We use the condition in Equation (5) to find which systems formed a CBD after the CEE phase. Motivated by the hydrodynamical simulations of Lau et al. [80] and the analysis of Tuna & Metzger [143], we take the fraction of envelope mass that forms the CBD (equation (1)) to be $f_{\mathrm{M}} = 0.05$, and adopt the specific angular momentum retention parameter (Equation (3)) to be $\beta= 3$ in our ‘fiducial CBD model’ (model 10). This also coincides with the observations of AGB disks that point towards high values of $\beta$ (Izzard & Jermyn [61]). Due to the large uncertainty in the amount of mass that remains bound after CEE and is available for CBD formation, we perform a parameter space exploration of $f_{\mathrm{M}}$, varying its value from 0.005 to 0.1 (models 11–22). We set its upper limit to keep the CBD gravitationally stable (e.g., Teasdale & Stamatellos [135]; see discussion in § 7.2.1). We also run additional models varying the value of $\beta$ from half to six (models 23–46), keeping $f_{\mathrm{J}} = \beta f_{\mathrm{M}}$ below one, as required by angular-momentum conservation. In our ‘fiducial CBD model’, we account for super-Eddington accretion using the prescription of Ghodla & Eldridge [42] (Equation (7)), taking $f_{\mathrm{GE}} = 2/3$, as inferred for NS accretors. To explore the sensitivity of our results to this treatment, we perform an additional simulation using the prescription of Tuna & Metzger [143] (Equation (8); model 47).

CBD Interactions in DNS Progenitors

Impact of CBD Mass and Angular Momentum

We study the interaction between CBDs formed after the CEE phase and the inner binaries in DNS progenitors. In the left panel of Figure 3, we show how the average fraction of the envelope accreted by the binary and its components varies with the fraction of the envelope that is redistributed into a CBD. We note that even when accounting for the suppression in the accretion rate due to super-Eddington outflows, the NS, which is the less massive component of the binary, tends to accrete a larger fraction of the CBD mass.

Let us explore our ‘fiducial CBD’ model (model 10), where 5% of the envelope remains bound and forms a CBD around the NS and the core, as an example. In this case, the central binary accretes on average 4% of the envelope, which is about 76% of the CBD mass. Taking an envelope mass of $\simeq10M_{\odot}$, motivated by the peak of its mass distribution (right panel of Figure 3), implies that about $0.38M_{\odot}$ is accreted by the binary system. Roughly $0.26M_{\odot}$ of this mass is accreted by the NS at a rate of $\simeq80$ times the Eddington rate, and $0.12M_{\odot}$ is re-accreted by the giant’s core that would later explode and form a DNS. This mass gain leads to heavier NSs, and, consequently, higher binary chirp masses (see § 6). The super-Eddington accretion rates onto the NS also raise the possibility that some of these post-CEE binaries could appear as ultraluminous X-ray sources during their interaction with the CBD (e.g., Kaaret et al. [64] for a review). However, their X-ray luminosities and observability would depend on the accretion flow geometry, beaming, and obscuration by surrounding material. Adopting the Tuna & Metzger [143] super-Eddington accretion allows for an average accretion rate by the NS of only five times the Eddington limit, reducing the amount of mass accreted by the NS to $\sim0.02M_{\odot}$, while keeping the mass accreted by the core mostly unchanged.

Besides driving gas accretion in which the binary components gain mass, the CBD influences the evolution of the binary through gravitational torques that exchange angular momentum with the orbit. In Figure 4, we show the CBD-induced change in the semimajor axis of DNS merger progenitors, combining all binaries that evolve into DNSs that merge within 14 Gyr either with or without CBD interactions in our fiducial models. As shown in the left panel, 90% of the DNS merger progenitors form CBDs after CEE, and most of them contract significantly as angular momentum is transferred from the binary to the disk. This reduces the number of systems that ultimately form DNSs, especially for tighter bina-

ries, as can be seen in the right panel. Although CBD interactions also significantly alter the post-CEE eccentricity, producing a mean increase of 0.16, this change is largely erased in our models by circularization after a case BB mass transfer episode that occurs between the core and the NS once the CBD is depleted.

Populations of Merging DNSs

We explore how the CBD-driven changes in the mass and orbital separation of DNS merger progenitors shown in Figure 3 and Figure 4 influence their subsequent evolution and final fate. As shown in Figure 5, CBD interactions prevent $\simeq29\%$ of the binaries that would otherwise have evolved to merging DNSs from doing so (left panel), while enabling the formation of new merging DNSs in systems that would otherwise result in different evolutionary outcomes, accounting for $\simeq23\%$ of the new DNS merger population (right panel).

In the presence of a CBD, many binaries that would have formed merging DNSs follow different evolutionary pathways instead (left panel of Figure 5). First and most importantly, the orbital shrinkage in some DNS progenitors leads to ‘premature’ mergers between the NSs and stellar cores mostly during the interaction with the CBD, preventing the binary systems from evolving to DNSs. These ‘premature’ mergers account for $\simeq81\%$ of the DNS mergers that no longer occur when CBD interactions are included, and can potentially power LFBOTs (see § 4). Another $\simeq9\%$ of the missing DNS mergers are instead lost to carbon–oxygen white dwarf (WD) – NS binaries. In these systems, CBD-induced orbital shrinkage initiates or enhances subsequent mass transfer from the stripped core to the first-born NS, removing sufficient mass from the core such that it forms a WD rather than undergoing the stripped-envelope SN that would otherwise produce the second NS. A quarter of these systems pass through an ultracompact helium-star-donor X-ray phase (e.g., [157]), before forming detached NS–carbon–oxygen WD binaries. None of the resulting WDs subsequently transfers mass onto the NS in our simulations. If these systems were later driven into contact, their relatively high WD masses would make stable, long-lived NS–WD ultracompact X-ray binary evolution unlikely ([11]). A remaining $\simeq6\%$ of binaries that would have evolved to DNS mergers without a CBD get disrupted during the second SN. In these systems, accretion from the CBD leads to secondaries with larger CO core masses (left panel of Figure 3) that receive larger natal kicks when exploding as stripped envelope SNe in our model, and are therefore more likely to disrupt binaries.

On the other hand, the loss of these DNS mergers is largely compensated by systems that form merging DNSs only when CBD interactions are included (right panel of Figure 5). In $\simeq44\%$ of these systems, CBD-driven orbital contraction leaves the binaries sufficiently tightly bound to survive the second SN and form DNSs in short-enough periods for future merger. Another $\simeq45\%$ originate from binaries that would otherwise have evolved to non-merging DNSs; their reduced separations shorten the GW inspiral time sufficiently for them to coalesce within the age of the Universe. In a $\simeq4\%$ minority, the CBD-induced orbital shrinkage leads to earlier mass transfer episodes between the stellar core and the NS in NS–BH progenitors, reducing the mass of the core and leading to the formation of an NS rather than a BH.

DNS Merger Rates and Timescales

To estimate volumetric event rates from our population synthesis data, we normalize the fraction of simulated binaries producing each event to the local star-formation-rate density, giving

where $f_{\mathrm{event}} = N_{\mathrm{event}}/N_{\mathrm{total}}$ is the number of events of interest divided by the total number of simulated binaries, $n_{\mathrm{SFR}} \simeq10^{7} M_{\odot}\,\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ is the local star formation rate normalized to the Kroupa IMF ([85]), and $M_{\mathrm{SF,sim}} \simeq88M_{\odot}$ is the total stellar mass formed per simulated massive binary. Substituting these values in Equation 9 results in

We present both the volumetric DNS merger rate $R_{\mathrm{DNS}}$ (left/brown $y$ axis) for the ‘fiducial CBD’ model, computed by substituting the simulated DNS merger fraction into Equation 10, and the relative reduction in the DNS mergers compared to the ‘fiducial no CBD’ model (right/blue $y$ axis) in Figure 6. For our ‘fiducial CBD’ case, where the specific angular momentum budget parameter $\beta$ equals 3 and CBDs that contain 5% of the envelope mass exchange angular momentum with the central binary, the interaction of the binary with the CBD leads to a $\simeq8\%$ reduction in the local DNS merger rate compared to a case where the formation of CBDs is not taken into account (from $R_{\mathrm{DNS}} = 23$ to $21\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$). The net rate change reflects the variety of competing evolutionary effects that reshape the population of merging DNSs, as illustrated in Figure 5. Overall, we find that CBD interactions prevent the formation

of more merging DNSs than they create, resulting in a lower merger rate.

Furthermore, we show the effect of varying different CBD parameters on the merger rate of DNSs. The left panel of Figure 6 presents its dependence on the CBD specific angular momentum budget parameter $\beta$. For most values of $\beta$, we find a decrease in the DNS merger rate relative to the no CBD case, driven by CBD-induced orbital shrinkage. This reduction is strongest for $\beta\lesssim1.9$. The reason is that, in this regime, fewer DNS progenitors form CBDs, and those that do are generally in tighter systems, making a ‘premature’ merger between the NS and the giant’s core more likely. For larger values of $\beta$, wider binaries that would otherwise evolve into non-merging DNSs form CBDs as well (see Unger et al. [145]), and the orbital shrinkage they induce can drive the DNS systems to merger. As a result, the number of mergers increases relative to the minimum, but still remains lower than in the no CBD case. When exploring how the disk mass affects the DNS merger rate (right panel of Figure 6), we find that increasing it generally lowers the rate density, as more massive disks exert stronger gravitational torques, driving greater orbital shrinkage (Figure 4) that makes the systems more prone to ‘premature’ mergers.

Under the Tuna & Metzger [143] super-Eddington prescription (model 47; Equation (8)), the population of merging DNSs changes in a manner similar to that shown in Figure 5. However, the overall DNS merger rate decreases only slightly, because the formation of new merging DNSs is approximately balanced by the loss of existing ones.

The interaction between the CBD and the inner binary affects not only the origin of the merging population and the total merger rate, but also the merger timescales of the systems that still manage to merge. Figure 7 presents the delay time distribution of DNS mergers for the fiducial models. We find that CBDs significantly reduce the number of DNS mergers prior to $\simeq0.2$ Gyr, lowering their merger rates below $\simeq40$ Myr by 40%. This deficit arises since the tightest post-CEE binaries, which would otherwise form DNSs with the shortest delay times, are the most likely to undergo ‘premature’ mergers between the NS and the companion’s core. This could influence r-process nucleosynthesis at early times, as discussed in § 7.1.

Based on these results, we conclude that CBDs tend to decrease the numbers of merging DNSs, especially at early times. Across the CBD parameter variations shown in Figure 6, the simulated rates we find span 17–23 Gpc$^{-3}$ yr$^{-1}$, well within the range inferred by The LIGO Scientific Collaboration et al. (2026) of 5.1–154.7 Gpc$^{-3}$ yr$^{-1}$, so current observations do not yet allow to discriminate different CBD parameters. We note, however, that the DNS merger rate inferred from the Galactic DNS population is somewhat higher; converting the Galactic estimate of Grunthal et al. [52] to volumetric units using a Milky Way star-formation rate of 1.65$M_\odot$ yr$^{-1}$ (Licquia & Newman [82]) and the same local star-formation-rate density adopted here gives $\simeq$ 140–310 Gpc$^{-3}$ yr$^{-1}$. Therefore, our absolute merger rates in both CBD and no-CBD models are lower than those favored by Galactic DNS observations, although these Galactic rates remain in some tension with the current GW-inferred rate.

CBD-induced LFBOTs

As discussed in § 3, the decrease in the DNS merger rate is primarily driven by ‘premature’ mergers between the NS and the giant’s core, caused by CBD-induced orbital hardening. Such mergers have been proposed as a possible channel for producing AT2018cow-like LFBOTs (Soker et al. [127]; Metzger [92]; Soker (2022); Tuna & Metzger [143]; Cohen & Soker [19]). LFBOTs are characterized by relatively short optical rise times of a few days and luminosities comparable to SN explosions (e.g., Margutti et al. [90]; Coppejans et al. [21]; Ho et al. (2020)), with evidence for dense circumstellar material that can include hydrogen-depleted gas (Fox & Smith [33]). In CBD-induced NS–core mergers, rapid accretion of the disrupted core onto the NS can power luminous emission through accretion-driven outflows, and the circumstellar material can consist of hydrogen-depleted gas lost by the stripped core at the onset of the merger, as well as hydrogen-rich material from the CBD and its photoevaporative winds [92, 143]. CBD-fed accretion onto the NS could also lead to the launching of jets that interact with the previously ejected envelope on a much shorter timescale, providing another route to powering LFBOTs (Soker 2022).

In Figure 8 we show the rate of CBD-induced NS–core mergers that can power LFBOTs as a function of our CBD parameters, computed by substituting their simulated number into Equation (10). We find that the rate of these ‘premature’ mergers initially rises with the specific angular momentum budget parameter $\beta$, before plateauing for $\beta\gtrsim2$ (left panel). At this point, most of the DNS progenitors that are close enough to experience an early merger do so, and after that, the wider systems in which CBDs are formed shift mostly from non-merging to merging DNSs. In contrast, the rate rises consistently with increasing CBD mass (right panel), as more massive disks exert stronger forces on the binary, and therefore drive more efficient orbital shrinkage. Moreover, we find that in the fiducial case about 87% of the ‘premature’ mergers happen within the CBD lifetime, allowing for a denser circumstellar material around the transient. Performing a similar analysis for binary BH progenitors, we find they have a negligible contribution to the population of LFBOTs.

Ho et al. (2023) compute the observed volumetric rate of LFBOTs based on data from the Census of the Local Universe survey and find it to be $R_{\mathrm{LFBOTs,obs}} \simeq70\ \mathrm{Gpc}^{-3}\ \mathrm{yr}^{-1}$, $\simeq15\%$ of which could be explained by CBD-induced NS–core mergers with fiducial model parameters. However, we note that the observed rate is highly uncertain, as combining estimates from Census of the Local Universe survey and the Bright Transient Survey give a conservative range of $R_{\mathrm{LFBOTs,obs}} \simeq0.3 - 420\ \mathrm{Gpc}^{-3}\ \mathrm{yr}^{-1}$ [21] (Ho et al. 2023). Therefore, a larger observational sample will be essential for assessing the contribution of CBD-induced NS–core mergers to the LFBOT population, which Rubin/LSST may provide by detecting tens to hundreds of LFBOTs per year [98]. Moreover, making a cleaner comparison between the simulated and observed rates would require integrating over the metallicity distribution representative of the local star-forming Universe.

In Figure 9 we present the post-first-SN systemic velocities of the LFBOT progenitors, $v_{\mathrm{sys}}$, vs the time elapsed between the first SN and the ‘premature’ merger of the NS with the stripped stellar core, $t_{\mathrm{d}} \equiv t_{\mathrm{LFBOT}} - t_{\mathrm{first\ SN}}$. Together, these quantities provide a first-order estimate of the typical offset of LFBOTs from their local stellar population following the first SN, due to the post-SN systemic velocity, $d \simeq v_{\mathrm{sys}}t_{\mathrm{d}}$, without taking the galactic potential into consideration. The combination of relatively small systemic velocities, $v_{\mathrm{sys}} \simeq2 - 27\ \mathrm{km\ s}^{-1}$, and short post-SN delay times, $t_{\mathrm{d}} \simeq0.01 - 7\ \mathrm{Myr}$, results in typically small displacements before the LFBOT occurs, with 40% having offsets smaller than 10 pc, and 99% smaller than 100 pc relative to the local stellar population. Therefore, we find that natal kicks are expected to make only a minor contribution to the observed locations of CBD-induced LFBOTs, and they are therefore expected to remain spatially associated with their local progenitor populations.

Furthermore, our results indicate that LFBOTs should occur near young or intermediate-age stellar populations, but not necessarily within regions of ongoing star formation, despite their small post-SN displacements. While massive stars are born in star forming regions, over the several tens of Myr preceding the LFBOT (color bar in Figure 9), the original H II regions are expected to have faded (e.g., [53]) and their natal molecular cloud to have been dispersed by stellar feedback (e.g., [16, 67]). This is consistent with the results of Nugent et al. [98], who find that although LFBOTs occur in actively star-forming galaxies, their locations do not always coincide with the brightest regions of their hosts. Events located several kilo-parsecs from any young stellar population would, however, be difficult to explain through the small post-SN displacements estimated from our model, and would require substantial pre-SN motion of the progenitor binary, potentially following dynamical ejection from its birth cluster (e.g., Oh & Kroupa 2016).

We emphasize that the displacements estimated here represent the offset produced following the first SN by the post-SN systemic velocity relative to the local stellar population, and are distinct from the offsets from host galaxy centers commonly reported for LFBOTs. Therefore, any large observed galactic offsets, such as $\simeq1.7\ \mathrm{kpc}$ for AT2018cow [106] and $\simeq16.5\ \mathrm{kpc}$ for AT2023fhn [17], primarily reflect the location of the progenitor population within its host galaxy. Moreover, pre-SN dynamical interactions, which are not included in our isolated binary evolution simulations, may also contribute to the observed displacements.

Sensitivity to additional population synthesis assumptions

Stellar Variations with CBD Interactions

In § 3.3, we have shown that accounting for the presence of a CBD in DNS progenitors can decrease the expected merger rate. Here, we test this result under additional variations of uncertain population synthesis parameters.

We evaluate the sensitivity of our findings to the underlying assumptions by varying model parameters related to the initial metallicity distribution (model 48), SN explosion prescription (models 49), mass transfer stability criterion (models 50-52), and CEE treatment (models 53-55), relative to the ‘fiducial CBD’ case, as described in § 2.2. We explore the effect of these variations on the DNS merger rate and the occurrence of LFBOTs, and find a nearly robust overall trend: CBDs lower the DNS merger rate (green bars in Figure 10), often by an amount comparable to or larger than the difference between our fiducial models. This reduction is mainly driven by ‘premature’, LFBOT-producing mergers (purple bars in Figure 10).

We note that in the Fryer et al. [35] case (models 3, 49), however, CBD interactions increase the rate of DNS mergers. The reason is that in this case CBDs produce two large competing effects. While they drive tight NS+stripped core binaries to ‘premature’ stellar merger before the second SN, as seen in the purple bars, they also create a larger amount of new DNS mergers. These mostly originate from previously non-merging and unbound DNS systems.

On the other hand, when adopting a CEE efficiency parameter of $\alpha= 2$ (model 8, 54), the decrease in the DNS merger rate due to CBD interactions is notably larger than in our ‘fiducial CBD model’. As in the fiducial models, in this case most of the missing DNSs are due to CBD-induced ‘premature’ mergers that lead to LFBOTs. Since for $\alpha= 2$ many more systems survive the CEE but end up as tight binaries (see § 5.2 for more details), the CBD-induced orbital shrinkage can drive the systems to an early merger more easily. As a result, the merger rate is reduced more strongly relative to the corresponding no-CBD model.

Stellar Variations without CBD Interactions

To assess the significance of our result given the current population synthesis accuracy, and to explore whether variations in the stellar-evolution prescriptions can affect the DNS merger rate in a similar way to CBDs, we perform simulations with the aforementioned variations relative to the ‘fiducial no CBD’ case (models 2-9). We find that changes in the assumed metallicity distribution, critical mass ratio for unstable mass transfer and post-CEE eccentricity have little to moderate effect on the DNS merger rate, whereas the treatment of SN prescription and common envelope removal can significantly affect the DNS merger rate, as shown in Figure 11.

We find that assuming a log-uniform metallicity distribution lowers the DNS merger rate by $\simeq17\%$ compared

to the ‘fiducial no CBD case’ in which all binaries are formed at solar metallicity, a result that is in line with recent studies finding a weaker metallicity dependence for DNS mergers than for BH-containing binaries (e.g., [43, 8, 130, 148]). Lower metallicity stars are typically more compact due to their lower opacity. The smaller radii at low metallicities tend to shift the onset of mass transfer to later evolutionary stages; as a result, stripping is less efficient, or can occur after the donor core has had more time to grow, implying larger core masses and a shorter Wolf–Rayet phase (e.g., [111, 68]). This can reduce the amount of wind-driven mass loss, which is also weaker at lower metallicities (e.g., [156, 155, 94]), and drive the formation of heavier SN progenitors. These higher-mass progenitors typically receive larger natal kicks within the Mandel & Müller [87] formalism, leading more binaries to unbind during the first SN, and therefore to the formation of fewer DNS systems, and fewer mergers as a result. Heavier progenitors also form more BHs. Another effect of smaller stars is that their envelopes are more tightly bound, giving an increase in stellar mergers during CEE. We note, however, that the effect of metallicity on the DNS merger rate is non-monotonic, with the largest reduction in rate occurring around $0.1Z_\odot$. Adopting the Fryer et al. [35] compact remnant mass function, together with a bimodal Maxwellian distribution for the imparted natal kick, increases the DNS merger rate by $\simeq90\%$. This is a result of a large reshuffling of the DNS merger population. In this model, the second SN that forms the DNS receives substantially lower natal kicks, leading to fewer binary disruptions at this stage, and consequently more systems that evolve to DNS binaries and later merge. In addition, the Mandel & Müller [87] prescription allows for BH formation in lower mass progenitors, turning some NS–BHs into DNS mergers.

We also explore how varying the CEE efficiency parameter $\alpha$ affects the DNS merger rate, and find that relative to the fiducial model with $\alpha= 1$, $\alpha= 2$ leads to a substantial increase in the DNS merger rate by about 65%. In this case, more energy is available to unbind the envelope, making envelope ejection easier and leading to a higher CEE survival rate, and consequently more DNS formation and mergers. This effect also allows heavier progenitors to form double compact objects, instead of ending their evolution as compact object-core mergers (e.g., Fryer & Woosley [36]; Soker & Gilkis [126]; Grichener & Soker [51]; Schrøder et al. (2020); Grichener [47]; Grichener [48]).

Shift in the Chirp-Mass Distribution of Merging DNSs

As in the case of binary BHs, the overall merger rate of DNSs is not sufficient to distinguish the underlying physics. In § 5.2 we found that varying other population synthesis assumptions, such as the SN prescription and the mass transfer stability criteria, could produce rate changes comparable to the inclusion of CBDs at the end of the CEE phase. Therefore, to provide an additional observational constraint, we turn now to the chirp masses of merging DNSs.

In the left panel of Figure 12 we show the pre-merger chirp-mass distribution of DNSs in our fiducial models. As discussed in § 3.1, the post-CEE core can re-accrete part of the envelope through a CBD. As a result, it can gain mass prior to exploding in the second SN event, usually leading to a heavier NS (see discussion on potential post-CBD mass loss in § 7.2.2). The first born NS can accrete mass as well, retaining material that is not launched through super-Eddington outflows. Therefore, when taking CBD interactions into consideration, the chirp-mass distribution shifts to larger values, with the prominent peaks moving by $\simeq0.05 - 0.1M_\odot$, and the mean value changing from $\simeq1.2M_\odot$ to $\simeq1.28M_\odot$. We find a deficit of mergers with $\mathcal{M} \lesssim1.08M_\odot$ compared to the models with no CBD, and a substantial excess of mergers with $\mathcal{M} \gtrsim1.35M_\odot$.

We next examine whether this shift toward higher chirp masses is specific to CBD interactions. For comparison, we consider the three additional non-CBD models that yield DNS merger rates comparable to that of the ‘fiducial CBD’ model (See Figure 11): log-uniform metallicity distribution (model 2), critical mass ratio for stable mass transfer from Claeys et al. [18] (model 4), and circular post-CEE eccentricity (model 9). The right panel of Figure 12 shows the resulting chirp-mass distribution of these models, and that of a CEE efficiency parameter $\alpha= 2$ (model 8). Although adopting the latter distribution yields a substantially higher DNS merger rate than the ‘fiducial CBD’ model, it is also the only non-CBD variation in our grid that results in a higher mean chirp mass ($\simeq1.24M_\odot$), which might be harder to distinguish from the effect of CBD-driven accretion. This arises since the higher CEE efficiency allows more massive progenitors to survive CEE and form DNSs, instead of merging within the envelope.

Therefore, we examine the chirp-mass of this case as well. We find substantial differences between the ‘fiducial CBD’ case and all presented models, which exhibit lower chirp-mass peaks, and significantly fewer systems with $\mathcal{M} \gtrsim1.35M_\odot$ than in the ‘fiducial CBD’ model.

To assess whether the structure in the predicted chirp-mass distributions is affected by the finite number of merging DNSs in our COMPAS populations, we bootstrap each model by resampling its merging-DNS population with replacement at fixed sample size and recomputing the chirp-mass PDF using the same bins as in Figure 12. The resulting 16th − 84th percentile intervals quantify the finite-sample uncertainty in each bin and are shown as error bars for the fiducial models in the left panel. We find that the broad features of the distributions and the shift toward higher chirp masses in models including CBD interactions are robust to these sampling uncertainties, while some of the smaller-scale variations between adjacent bins may arise from finite-sample fluctuations.

Unfortunately, only two DNS mergers have been detected so far, GW170817 with $\mathcal{M} = 1.188^{+0.004}_{-0.002}M_\odot$ ([1]) and GW190425 with $\mathcal{M} = 1.44^{+0.02}_{-0.02}M_\odot$ (Abbott et al. 2020), so it is not yet possible to construct an observational chirp-mass distribution for a meaningful comparison with our models. We note, however, that the high chirp mass of GW190425 is much more compatible with models that include CBD interactions. As in previous population-synthesis studies (e.g., Osłowski et al. [102]; Vigna-Gómez et al. [153]; Mandel et al. [88]; Chattaraj et al. [13]), this mass remains challenging to explain within most of our no CBD models (model 1-7, 9), for which fewer than 1% of merging DNSs have chirp masses in the range $1.42M_\odot\le\mathcal{M} \le1.46M_\odot$. However, when CBD-driven accretion is included, the fraction of DNS mergers with GW190425-like chirp masses increases to $\simeq10\%$, making such systems considerably more common in the simulated population. By contrast, the $\alpha= 2$ model (model 8), despite having the second-highest mean chirp mass, produces GW190425-like DNS mergers at a fraction lower by a factor of two relative to the ‘fiducial CBD model’. We note, however, that allowing super-Eddington accretion during the CEE phase or the subsequent case BB mass transfer could result in higher chirp masses as well.

For additional observational context, we also show in Figure 12 the eight Galactic DNS systems with measured individual masses that are expected to merge within a Hubble time (see Appendix § 12 for more details). The chirp masses estimated from the component mass measurements are in the range $\mathcal{M} \simeq1.10 - 1.25M_\odot$. While four Galactic DNSs lie in regions populated comparably by the fiducial CBD and no-CBD models, for the other four the no-CBD model predicts a higher local probability density. To date, no merging Galactic DNSs have been observed near the high-mass peak of the ‘fiducial CBD’ chirp-mass distribution. We note, however, that more massive Galactic DNSs tend to have shorter GW inspiral times (Grichener et al. [49]), and are therefore less likely to be represented in the present-day Galactic DNS population.

The substantial mass accretion onto NSs in our ‘fiducial CBD’ model (§ 3.1), which might explain the high mass of GW190425, is difficult to reconcile with the measured component masses of Galactic DNSs, suggesting that CBDs may form in a smaller fraction of systems and/or that less mass is accreted from them than found in our ‘fiducial CBD’ model. However, the ‘fiducial no CBD’ model also fails to reproduce the observed component masses, indicating that neither model fully captures their distribution. Interestingly, examining the individual masses reveals that in most systems, the recycled pulsar, generally identified as the first-born NS, is the more massive component. This could potentially be explained by additional mass gained through accretion from a CBD.

Discussion

R-process Enrichment

It remains unclear whether DNS mergers alone can account for all r-process nucleosynthesis (e.g., Ji et al. 2019; Holmbeck & Andrews [58]; Soker 2026). Recent studies have suggested that due to their long GW inspiral times, DNS mergers might not occur early enough to explain the production of heavy elements in the early Galaxy (e.g., Côté et al. 2019; Wehmeyer et al. [159]; Kobayashi et al. [69]; Kobayashi et al. [70]; Chen et al. [14]; Saleem et al. [114]). The strong decrease in the DNS merger rates at timescales of $\lesssim40$ Myr due to CBD-induced ‘premature’ mergers (Figure 7) makes this tension even more severe. Such a substantial suppression of short-delay mergers strengthens the need for an additional, more prompt r-process source, and may extend the time over which such a source is needed.

On this timescale, which is comparable to the lifetime of massive stars, r-process nucleosynthesis could be dominated by the merger of an NS with the core of a giant star during CEE (e.g., Papish et al. [104]; Grichener & Soker [51]; Grichener et al. [50]; Grichener [47]; Jin & Soker [63]; Grichener [48]; Soker 2025; Soker 2026; Thielemann & Cottam [138]). In this case, the r-process nucleosynthesis occurs in the self-neutronizing accretion disk that is formed around the NS from the destroyed core material. These mergers are significantly faster than DNS mergers since they are dominated by the dynamical friction of the envelope rather than caused by GW emission. As CBD interactions drive a large fraction of DNS progenitors to undergo a ‘premature’ merger between the NS and the stellar core, thereby introducing an additional NS–core merger channel, their relative contribution to r-process nucleosynthesis may become even more significant, as will be explored in a future study.

Limitations of the Physical Models

Disk Properties Uncertainties

Mass: The amount of envelope mass that remains bound after CEE and can subsequently be redistributed into a CBD varies significantly across hydrodynamical simulations, with some studies suggesting that the majority of the envelope is ejected (e.g., Lau et al. [80]), and others indicating most of it may remain bound (e.g., Ohlmann et al. [100]). Since the effects of the disk self-gravity are not included in the hydrodynamical simulations on which our model is based, when varying the CBD mass in our population-synthesis study, we take 10% of the envelope mass as an upper limit. For this value, the Toomre parameter (e.g., Toomre [140]; Kratter & Lodato [72]) satisfies $Q \gtrsim1$ when assuming a constant aspect ratio of $h = 0.1$, and the disks are expected to be marginally stable against local gravitational instability. Higher CBD masses would reduce $Q$, increasing the importance of self-gravity, and potentially leading to clumping or fragmentation if the disk cools sufficiently efficiently (e.g., [39, 73, 112]). However, due to uncertain disk properties, including its temperature, surface-density profile, and radial extent, effects of self-gravity cannot be ruled out for all systems within our explored parameter space.

Lifetime: The CBD simulations performed by Siwek et al. [120] provide the mass accreted towards the central binary, but are agnostic to the accretion rate and the disk lifetime, as the effect of photoevaporation and disk wind are not included. We take a characteristic lifetime of $\tau_{\mathrm{CBD}} = 10^{5}\ \mathrm{yr}$, comparable to the thermal timescale of a stripped helium core and, in some cases, even to the nuclear timescale of a more evolved core. For a fixed CBD mass, a shorter disk lifetime implies a higher initial accretion rate. In most of the models explored here, which adopt the super-Eddington correction of Ghodla & Eldridge [42], the accretion timescale does not affect the total mass accreted once the flow becomes super-Eddington (Equation (7)). By contrast, under the prescription of Tuna & Metzger [143], higher accretion rates drive stronger super-Eddington winds and may therefore result in less mass being accreted onto the NS (Equation (8)). The resulting NS mass growth can therefore depend appreciably on the assumed disk lifetime. However, because the total NS mass gain is already small under this prescription, this additional variation has little impact on the population-level outcomes considered here. Reducing the disk lifetime by an order of magnitude, from $\tau_{\mathrm{CBD}} = 10^{5}\ \mathrm{yr}$ to $\tau_{\mathrm{CBD}} = 10^{4}\ \mathrm{yr}$, changes the merger rates explored here by less than 3%.

Thickness and Viscosity: There are additional disk properties that might affect the interaction of a CBD with a binary system. The aspect ratio (set to $h = 0.1$ for the CBDs explored in Siwek et al. [120]) and viscosity (treated through the disk prescription of Shakura & Sunyaev [118] with $\alpha_{\mathrm{visc}} = 0.1$) of the disk play a crucial role in shaping its interaction with the inner binary and its subsequent orbital evolution, as they control the structure of the inner cavity and the accretion streams. Thicker disks have greater pressure support that smooths the cavity edge, weakening the tendency toward orbital contraction and, in many cases, instead driving expansion of the binary. Thinner disks preserve a sharper, more strongly non-axisymmetric cavity morphology, which enhances the gravitational pull from disk gas near the binary and can drive stronger orbital contraction (e.g., [139, 27, 28]). The viscosity then modulates how efficiently gas is transported across the cavity edge: higher viscosity refills and mixes the inner region more effectively, lowering density contrasts and weakening the gas-driven force, whereas lower viscosity maintains a deeper, higher-contrast cavity with more coherent, lagging gas streams, strengthening the gravitational torque that removes angular momentum from the binary and driving stronger orbital decay. Moreover, the adopted viscosity prescription might influence the effect of CBDs on the eccentricity evolution of the binary in our estimated post-CEE eccentricity regime (see § 7.2.2). D’Orazio & Duffell [29] find that equal-mass binaries with initial eccentricities $e \lesssim0.1$ are circularized by the CBD, while the simulations of Siwek et al. [120] underlying our framework instead predict eccentricity growth for equal-mass binaries with $e \gtrsim0.01$. This disagreement illustrates that the eccentricity evolution, particularly for low eccentricity binaries, remains sensitive to the adopted disk physics. The thickness and viscosity of post-CEE disks are not observationally constrained and might vary throughout the disk lifetime, with the disk being initially thicker just after the CEE and becoming thinner as it cools. While the values in the simulations of Siwek et al. [120] are in the range of commonly used values for studying CBDs formed after CEE (e.g., [143, 160]), a parameter-space exploration of these values would be valuable for testing the robustness of our results, particularly as preparation for future observations that may help constrain their relevant range.

Stellar Evolution Uncertainties

Effect of Metallicity: In our fiducial model we assume that binaries form at solar metallicity, and relax this assumption by adopting a log-uniform distribution in metallicity between $10^{-4}$ and 0.03. In reality, DNS progenitors span a range of metallicities that depends on the star formation history and redshift of their host galaxy. Since the metallicity affects the stellar radii and wind mass loss rates (as discussed in § 5), it would influence the orbital evolution. Most cosmic star formation occurs at relatively high metallicities (e.g., [147]). Lower metallicities are found at higher redshifts, but unless there are large delays between formation and merger, mergers in these environments are generally harder to detect (e.g., [1]). Therefore, we adopt solar metallicity as a representative choice for the distributions most likely to be observed in the near future. A more complete treatment, involving simulations across a metallicity grid and convolution with a cosmic star formation and metallicity evolution history, including selection effects, is beyond the scope of this work.

Post-CEE Circularization: Post-CEE systems typically exhibit small eccentricities (e.g., [9, 75]), and population synthesis studies therefore commonly assume circular orbits after CEE. However, complete circularization during CEE is not guaranteed, and hydrodynamical simulations can leave binaries with non-zero residual eccentricities (e.g., [105, 44]). Since CBD interactions and especially eccentricity pumping can be highly sensitive to the initial eccentricity of the binary, this uncertainty is important for our models. Moreover, the CBD-driven evolution of binaries at low but non-zero eccentricities is itself uncertain. Within SPINDLER, binaries starting from exactly circular orbits do not develop eccentricity. We therefore adopt a simple phenomenological prescription in which after CEE the binaries retain 10% of their pre-RLOF eccentricity in most models. We note that pre-RLOF tidal evolution was not included in our calculations (see e.g., Kapil et al. [65] for a new prescription), so this assumption might overestimate the post-CEE eccentricity. Previous work, however, suggests that a significant fraction of DNS progenitors may retain non-negligible eccentricities at the onset of CEE, even when accounting for tidal circularization (e.g., [154]). Nevertheless, within our adopted prescription, the binaries develop non-negligible eccentricities in the presence of a CBD, while the DNS merger rate remains similar to that obtained in simulations that assume CBD interactions and $e_{\mathrm{post,CEE}} = 0$. For even more eccentric post-CEE orbits, the Siwek et al. [120] prescription tends to enhance the eccentricity further, driving it towards the $e \sim0.5$ attractor they found (see also Figure 1 in Valli et al. [146]), and thereby strengthening its impact on the merger rate and associated observables.

Effect of CBD Torques on Mass Transfer: While our model includes CBD-driven orbital evolution, we rely on the standard COMPAS prescriptions to determine mass-transfer stability, without explicitly accounting for the influence of ongoing disk torques. Continued angular momentum extraction from the binary by the disk could enhance the mass-transfer rate and destabilize mass transfer that would otherwise remain stable (e.g., [92, 143]). This may be particularly relevant for case BB mass transfer, which is assumed to be stable in our simulations, potentially increasing the rate of NS–core mergers and associated LFBOTs, while further reducing the formation of merging DNSs. Quantifying the impact of disk torques would require a more self-consistent treatment of mass transfer in the presence of a CBD.

Post-CBD Mass Transfer: In our models, stripped stars re-accrete hydrogen-rich envelope material through a CBD. This may cause them to re-expand, as in the case of stars retaining residual envelopes following CEE (e.g., [77, 152]). Additional expansion may be driven by gravitational energy release and thermal readjustment during accretion ([79, 20]). Consequently, these stars might fill their Roche lobes, limiting further accretion and potentially redirecting material towards the NS companion, therefore altering the component masses and orbital evolution. These effects are not included in our models, and are left for future work.

Summary

In this study, we developed a framework to model CBD formation in post-CEE systems and its impact on the evolution of the inner binaries, and implemented this framework in the population synthesis code COMPAS (§ 2). We used it to explore the effect of CBD interactions on the population of merging DNSs and LFBOTs (§ 3–§ 6), and discussed the implications for r-process enrichment and the main limitations and caveats of our approach (§ 7).

To determine which systems form a CBD from leftover envelope gas surrounding the first-born NS and the exposed core, we followed Unger et al. [145] and compared the specific angular momentum available for CBD formation (equation (3)) to the Keplerian specific angular momentum required to support the disk (equation (4)). If the former is larger than the latter, we assumed the bound material has enough angular momentum to form a centrifugally supported CBD around the binary (equation (5)). In this case, we used SPINDLER (Valli et al. [146]), which interpolates the hydrodynamical CBD simulations of Siwek et al. [120], to model the CBD-driven evolution of the semimajor axis, eccentricity, and component masses of the binary shaped by gravitational and accretion forces. For compact objects, we reduced the amount of mass that is ultimately accreted and accounted for the angular momentum carried away by super-Eddington winds (equations (7)–(8)).

We used COMPAS to generate populations of massive binaries and to quantify how CBD interactions affect DNS mergers. We ran a suite of models, varying metallicities, SN and mass transfer treatments, and CBD parameters (Table 1). Our main findings are as follows.

1. CBDs change the merging DNS population and lead to a decrease in the merger rate: The interactions between the CBD and the central binary

typically shrink the orbit (Figure 4), changing the binary population that results in DNS mergers (Figure 5). This occurs mostly through ‘premature’ mergers between the NS and the giant’s core in systems that would otherwise evolve into merging DNSs. The reduction in the number of DNS mergers is partly compensated by systems that would be disrupted during the second SN or form DNSs too wide to merge within the age of the Universe, but instead become merging DNSs due to CBD-driven orbital contraction. On balance, however, we found that CBD interactions reduce the predicted merger rate relative to analogous models without CBDs in most cases (Figure 6, Figure 10). In our fiducial model, the DNS merger rate decreases by 8%, while in the most strongly affected case it can be reduced by as much as 28%.

(ii) The missing DNSs can power LFBOTs: The ‘premature’ NS–core mergers, driven by orbital shrinkage during the binary-CBD interaction, generally occur before the disk material is depleted and can therefore give rise to LFBOTs (§ 4). This links DNS mergers and LFBOTs through a common progenitor pathway. We found that CBD-induced LFBOTs occur at a fiducial rate of $\simeq11\ \mathrm{Gpc}^{-3}\ \mathrm{yr}^{-1}$, and are robust across our population study (Figure 8, Figure 10). They are expected to remain spatially associated with their local progenitor environment, with natal kicks causing only a minor displacement relative to the local stellar population (Figure 9). This implies the transients typically arise near relatively young stellar populations, without necessarily coinciding with sites of ongoing star formation.

(iii) CBDs preferentially suppress DNS mergers at short delay times: Since most DNS progenitors that undergo CBD-induced NS–core mergers originate from initially tight binaries, CBD interactions mostly suppress DNS mergers at short delay times, lowering their rate by 40% for delay times below $\simeq40\ \mathrm{Myr}$ (Figure 7). This may further challenge the ability of DNS mergers to account for early r-process enrichment (§ 7.1). However, these ‘premature’ mergers provide a new avenue for forming NS–core mergers, which might contribute to r-process nucleosynthesis on timescales comparable to those of massive-star evolution.

(iv) CBD interactions produce more massive merging DNSs: CBDs influence DNS progenitors by mass accretion, which shifts the chirp-mass distribution of merging DNSs to higher values (Figure 12), with an average fiducial growth of $\simeq0.08M_{\odot}$. This may help explain the high mass of GW190425, which is not accounted for by standard population models to date (§ 6).

While the parameter space exploration showed that our results are robust to most stellar evolution variations, the uncertainties in CBD parameters could still significantly affect our predictions. The thickness and viscosity of the disk, for instance, could alter the amount of mass that is accreted by the central binary, and lead to a stronger inspiral. A natural follow-up, beyond the scope of this paper, would be to relax our assumptions about disk thickness, either by interpolating over a broader suite of hydrodynamical simulations that spans a wider range of aspect ratios (e.g., Tiede et al. [139]; Dittmann & Ryan [27]), or by adopting an analytic CBD framework in which these parameters can be varied directly (e.g., Izzard & Jermyn [61]; Tuna & Metzger [143]; Wei et al. [160]). Combined observations of DNS merger properties, LFBOTs, and r-process enrichment could provide an indication of the presence of CBDs in post-CEE massive binaries.

In this work we focused on the impact of CBDs formed after the CEE phase on the population of merging DNSs, since the majority of DNS progenitors are expected to undergo at least one common envelope phase (e.g., Ta uris et al. [132]; Vigna-Gómez et al. [153]; Grichener [47]; Nathaniel et al. [96]). More broadly, the same framework can be applied to any type of stellar binary in which a CBD forms following CEE, with potential implications for post-asymptotic giant branch systems, luminous red novae, interacting binaries that contain WDs, and binary BHs.

Acknowledgments

We thank Mike Lau, Michaela Mapelli, Smadar Naoz, Fabian Schneider, and Fritz Röpke for helpful discussions. A.G. acknowledges support from the Steward Observatory Fellowship in Theoretical and Computational Astrophysics, the IAU-Gruber Fellowship, and the CHE Fellowship. R.V. was supported by the NOMIS-Fellowship Program. M.R. acknowledges support from NASA (ATP: 95380NSC24K0932). I.M. acknowledges support from the Australian Research Council (ARC) Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE230100016.

Data Availability

The population synthesis data generated for this project is publicly available on Zenodo at 10.5281/zen-do.23851460.

Software

Simulations in this paper used a modified version of the open-source rapid binary population synthesis code COMPAS, incorporating the CBD framework developed in this work. It is publicly available on GitHub. The CBD framework is integrated in COMPAS, starting from version 03.31.00. This work also made use of SPINDLER [146].

Galactic merging double neutron star observations

Galactic DNSs that are expected to merge within a Hubble time provide an additional observational comparison to the chirp-mass distributions predicted by our population models. In Table A1, we list the properties of the eight observed merging Galactic DNSs shown in Figure 12. Their chirp masses are calculated from the measured component masses, while their remaining GW inspiral times are estimated from the present-day orbital parameters according to the Peters [107] equations.

Name	$P_{\mathrm{m}}$	$e_{\mathrm{m}}$	$M_{\mathrm{p}}$	$M_{\mathrm{c}}$	$\mathcal{M}$	$\tau_{\mathrm{GW}}$	Ref
	[days]		[$M_{\odot}$]	[$M_{\odot}$]	[$M_{\odot}$]	[Gyr]
J1946+2052	0.078	0.064	1.2842	1.2477	1.102	0.044	(1)
J1757–1854	0.183	0.606	1.3384	1.3945	1.189	0.075	(2)
J0737–3039A/B	0.102	0.088	1.338185	1.248868	1.125	0.085	(3)
B1913+16	0.323	0.617	1.438	1.390	1.231	0.301	(4)
J1913+1102	0.206	0.090	1.599	1.290	1.249	0.467	(5)
J0509+3801	0.380	0.586	1.36	1.45	1.222	0.575	(6)
J1756–2251	0.320	0.181	1.341	1.230	1.118	1.66	(7)
B1534+12	0.421	0.274	1.3330	1.3455	1.166	2.74	(8)

References: (1) Meng et al. 2025; (2) Singha et al. 2026; (3) Kramer et al. 2021; (4) Weisberg & Huang 2016; (5) Miao et al. 2026; (6) Lynch et al. 2018; (7) Ferdman et al. 2014; (8) Fonseca et al. 2014.