1. Introduction
White dwarfs (WDs) are the most abundant stellar remnants constituting a significant component of the stellar population in the Universe. This makes them crucial probes for studying the structure and evolution of compact objects (Fontaine, Brassard, & Bergeron Reference Fontaine, Brassard and Bergeron2001; Althaus et al. Reference Althaus, Córsico, Isern and García-Berro2010). A considerable fraction of WDs exist in binary systems (a system comprising of two WDs), commonly referred to as double WDs (DWDs), which play a pivotal role in several astrophysical contexts. These systems often undergo multiple distinct mass transfer phases, making them ideal laboratories for investigating binary stellar evolutions. Furthermore, the coalescence of DWDs within such systems is one of the leading progenitor models for Type Ia supernovae (Webbink Reference Webbink1984; Meng & Yang Reference Meng and Yang2010; Sato et al. Reference Sato, Nakasato, Tanikawa, Nomoto, Maeda and Hachisu2015; Canals, Torres, & Soker Reference Canals, Torres and Soker2018). Additionally, DWDs are expected to dominate the low-frequency gravitational wave (GW) spectrum, potentially generating a stochastic background that can obstruct other GW source signals (Evans, Iben, & Smarr Reference Evans1987; Nelemans, Yungelson, & Portegies Zwart Reference Nelemans, Yungelson and Portegies Zwart2001c; Lamberts et al. Reference Lamberts, Blunt, Littenberg, Garrison-Kimmel, Kupfer and Sanderson2019). Moreover, Galactic DWD binaries create a dominant GW foreground for Laser Interferometer Space Antenna (LISA) with mass-transferring systems significantly enhancing the noise at higher frequencies while still leaving thousands of resolvable sources (Ruiter et al. Reference Ruiter, Belczynski, Benacquista, Larson and Williams2010). Finally, analysing close DWDs, especially considering the binary evolution theory, provides a way to test and refine cooling models for low-mass WDs (Sarna et al. Reference Sarna, Ergma and Gersˇkevitsˇ-Antipova2000; Istrate et al. Reference Istrate, Marchant, Tauris, Langer, Stancliffe and Grassitelli2016).
The formation of DWDs has been studied extensively using both analytical (Iben & Tutukov Reference Iben1986, Reference Iben1987) and numerical population synthesis methods (Lipunov & Postnov Reference Lipunov and Postnov1988; Yungelson & Tutukov 1993; Han et al. Reference Han, Podsiadlowski and Eggleton1995; Iben, Tutukov, & Yungelson 1997; Han Reference Han1998). These studies primarily differ in their assumptions about key evolutionary processes, such as common envelope (CE) ejection, mass transfer stability, and angular momentum loss, which lead to variations in predicted population characteristics (Korol et al. Reference Korol, Hallakoun, Toonen and Karnesis2022b; Scherbak & Fuller Reference Scherbak and Fuller2023; Toubiana et al. Reference Toubiana, Karnesis, Lamberts and Miller2024; van Zeist et al. Reference van Zeist, van Roestel, Nelemans, Eldridge, Korol and Toonen2025). A cross-comparison of these models provides valuable insights into how such modelling choices affect the predicted properties of DWD populations.
In recent years, the catalogue of known DWDs has grown substantially owing to surveys such as the Extremely Low Mass survey and Double-Lined DWD survey, which have collectively identified DWD systems with well-constrained masses and orbital periods (Brown et al. Reference Brown2020; Munday et al. Reference Munday2024, Reference Munday2025). This expanded dataset allows for a more precise comparison between theoretical predictions and observed populations. In particular, the discovery of ultracompact systems, such as J0526+5934 with an orbital period of just 20.5 min (Rebassa-Mansergas et al. Reference Rebassa-Mansergas2024), has enabled stringent tests of binary evolution at the shortest orbital separations.
In parallel, Gaia-based analyses have improved estimates of the local DWD fraction, revealing that approximately 1.2% of resolved WDs within 100 pc reside in DWD systems, while approximately 6.3% are in resolved binaries with the companion being a main-sequence star (Torres et al. Reference Torres, Canals, Jime´nez-Esteban, Rebassa-Mansergas and Solano2022). Furthermore, forward modelling efforts for LISA GW mission predict that up to approximately 16 000 Galactic DWDs may be individually detectable over a 10 yr observing run (Korol et al. Reference Korol, Hallakoun, Toonen and Karnesis2022b), offering a unique opportunity to study the compact binary population through GW observations.
In this study, we utilise the binary population synthesis code compas (Compact Object Mergers: Population Astrophysics and Statistics)Footnote a to investigate the influence of different binary evolutionary scenarios on the formation properties and accretion histories of DWDs. By systematically exploring key physical processes, such as Roche lobe overflow (RLOF), mass transfer stability, and core composition of WDs, we assess their impacts on the resulting WD core types, orbital properties, and mass distributions. In particular, we recover the previously reported bimodal separation distribution (Korol et al. Reference Korol, Hallakoun, Toonen and Karnesis2022b,a) with intermediate under-density in DWD orbital separations, and we extend earlier work by (i) providing a unified analysis that links formation channel (RLOF-only vs. CE-involving), WD core type (He/CO/ONe), and mass transfer rate within a single compas-based framework, and (ii) quantifying how the separation gap responds to changes in the CE structural parameter. Our results provide a physically interpretable map between the observed properties of DWDs and their underlying interaction histories, and thus offer a concrete framework for interpreting future detections, particularly in the context of upcoming space-based GW observatories like LISA.
The remainder of this paper is organised as follows. In Section 2, we review the key aspects of binary stellar evolution and the formation of DWDs. Section 3 provides our key results, and Section 4 discusses the results and compares them with the previous literature. Finally, concluding remarks are provided in Section 5.
2. Revisiting binary evolution physics
The evolutionary pathways of binary stellar systems are similar to those of single stars when their orbital separations are sufficiently wide that neither star fills its Roche lobe during their evolution. In this case, the stars evolve independently, without undergoing mass transfer or tidal interactions. Therefore, in wider binaries, each star evolves independently according to the single-star evolution process. However, in close binaries (systems with separations small enough that at least one star is expected to fill its Roche lobe within a Hubble time), gravitational interactions make their evolutionary paths interdependent, leading to complex phenomena that significantly reshape the evolution of both stars and the system itself.
A fundamental process in close binary stellar evolution is mass transfer, which is initiated when one star expands to fill its Roche lobe, a region gravitationally bound to the star. The Roche lobe radius for the donor star (
$R_L$
) can be approximated by the following relation given by Eggleton (Reference Eggleton1983) as
\begin{equation} \frac{R_L}{a} = \frac{0.49q^{2/3}}{0.6q^{2/3} + \ln\!({1 + q^{1/3}})},\end{equation}
where a is the semi-major axis of the binary and
$q$
is the mass ratio of donor to accretor.
The stability of mass transfer depends critically on the structure of the donor star and mass ratio. Stable mass transfer occurs if the donor star and its companion have comparable masses or the donor has a radiative envelope. In this case, matter is transferred through the inner Lagrangian point and accreted by the companion over a long timescale. This scenario allows both stars to adjust with the changing mass ratio without significant disruption. On the other hand, unstable mass transfer occurs if the donor is a convective giant or significantly more massive than its companion. In this case, the companion cannot efficiently accrete the incoming material, leading to envelope engulfment (Paczynski Reference Paczynski, Eggleton, Mitton and Whelan1976) and CE phase. However, the stability of mass transfer remains one of the open problems in binary evolution. The location of the stability boundary is governed by multiple physical factors, such as the donor’s internal structure (whether it has a radiative or convective envelope, or is partially degenerate), the binary mass ratio, the timescale of mass transfer (nuclear, thermal, or dynamical), and a number of poorly constrained processes such as angular momentum loss via outflows, tides, or magnetic braking. Recent studies have shown that even moderate changes in these assumptions can substantially alter the predicted stability threshold, particularly for giants and partially stripped progenitors. Consequently, the standard guideline of that donors with radiative envelopes undergo stable mass transfer while those with convective envelopes experience unstable mass transfer is now recognised as an oversimplification (Ge & Han Reference Ge and Han2024). This uncertainty in the stability boundary implies that the relative rates of stable RLOF and CE evolution in any rapid population synthesis model should be interpreted as model-dependent rather than as precise predictions.
In this work, we adopt the default compas stability criteria without recalibration, which are based on the rapid binary evolution models by Hurley et al. (Reference Hurley, Tout and Pols2002). Our population-level predictions are therefore rely on these prescriptions. We partially address this model dependence by studying the sensitivity of our results to the changes in significant parameters such as the degree of non-conservative mass transfer and the angular momentum lost through outflows. These variations can alter the occurrence of stable RLOF versus CE evolution. A comprehensive treatment, which would involve the combination of complete stellar-structure models with the mass transfer process on thermal and dynamical timescales, is not within the realm of a rapid population synthesis study. Nevertheless, we discuss these outcomes which remain robust across a plausible range of input parameters. Therefore, our research should be interpreted as an exploratory study of how different interaction channels imprint themselves on DWD population, rather than as a fully calibrated model of the observed systems.
The change in orbital separation during the CE phase is typically modelled using the energy formalism (Webbink Reference Webbink1984), given by
\begin{equation} \alpha_{\text{CE}}\left(\frac{GM_{\text{c,1}}M_{2}}{2a_{\text{f}}} - \frac{GM_1M_2}{2a_{\text{i}}}\right) = \frac{GM_1M_{\text{env,1}}}{\lambda_{\text{CE}} R_1},\end{equation}
where
$\alpha_{\text{CE}}$
is the CE efficiency,
$M_{\text{c,1}}$
is the core mass of the donor,
$M_{\text{env,1}}$
is the envelope mass,
$a_{\text{i}}$
and
$a_{\text{f}}$
are the initial and final semi-major axis,
$R_1$
is the radius of the donor,
$M_1$
and
$M_2$
are the mass of the donor and the companion, and
$\lambda_{\text{CE}}$
is the parameter that describes the structure and binding energy of the donor star’s envelope. While the energy formalism is commonly applied in rapid population synthesis, alternate models exist. The angular momentum based prescription proposed by Nelemans et al. (Reference Nelemans, Verbunt, Yungelson and Portegies Zwart2000) links orbital contraction to angular momentum loss rather than energy considerations. These two formalisms represent competing physical interpretations of the CE process, and the choice between them can substantially alter the predicted post-CE orbital separation distribution. In the present work, we adopt only the energy formalism and do not implement the angular-momentum prescription, so our post-CE separation distributions should be interpreted within this framework. During this brief and dynamic CE phase, drag forces within the envelope cause the stellar cores to spiral inward, which may lead to ejection of the envelope and formation of a close binary. If the envelope ejection is inefficient, it may lead to a stellar merger. Thus, the outcome of this CE phase determines whether the system evolves into a compact binary or coalesces.
These complex evolutionary pathways lead to the formation of diverse compact binary systems and unusual stellar phenomena, including cataclysmic variables, X-ray binaries, Algol-type systems, and blue stragglers, which often cannot be explained by a single stellar evolution alone (Pringle & Wade Reference Pringle and Wade1985; Wijers, Davies, & Tout 1996; Perets & Fabrycky Reference Perets and Fabrycky2009. Although many details of close binary evolution are uncertain, a robust theoretical framework now exists that enables us to track how binaries progress through different evolutionary phases. In the context of DWD formation, progenitor properties, such as stellar mass, core composition, orbital period or semi-major axis, and eccentricity, determine the sequence of interactions that may eventually produce a binary system containing two WDs. Mapping of these evolutionary pathways is essential for interpreting the observed DWD population and predicting the outcomes of future GW surveys.
2.1. Double white dwarf formation
Close DWDs, which in our simulations correspond to orbital separations of
$\sim0.05-100\,\mathrm{R}_\odot$
(i.e. periods from minutes to days/months), are predominantly formed through two widely discussed families of evolutionary channel. In the first case, the primary star evolves to fill its Roche lobe, initiating stable mass transfer to the secondary star, which may or may not accrete the material. This phase produces an intermediate binary consisting of the first-formed WD and a secondary star (main-sequence or asymptotic giant branch). Subsequently, the secondary star fills its Roche lobe, triggering a second mass transfer episode. If this second mass transfer phase remains dynamically stable, the initially more massive primary has already become a WD, so mass is now transferred from the more massive secondary to the lighter WD (a mass ratio greater than unity). This configuration generally tends to shrink the orbit, and in our non-conservative models, additional angular momentum is carried away by the expelled material. As a result, binaries that undergo two stable RLOF episodes (RLOF+RLOF) can also produce close DWDs with separations of only a few solar radii, even without ever entering a CE phase. If the donor is a convective giant or significantly more massive than the WD companion, this phase becomes dynamically unstable (Paczynski Reference Paczynski, Eggleton, Mitton and Whelan1976; Tutukov & Yungelson Reference Iben1997; Webbink Reference Webbink1984; Iben & Livio Reference Iben1993). During this phase, the core of the donor and its companion spiral together within the shared envelope because of drag forces, which release orbital energy. This energy heats and ultimately expels the envelope, resulting in a significant reduction in the orbital separation. This explains how DWDs with separations of only a few solar radii form from initially wider binaries (Woods et al. Reference Woods, Ivanova, van der Sluys and Chaichenets2012; Ivanova et al. Reference Ivanova2013).
In the other scenario, both stars undergo consecutive CE phases, each initiated by an unstable RLOF. In our rapid binary-evolution framework, every CE event is triggered by dynamically unstable mass transfer (or an equivalent envelope-engulfment instability), so episodes of unstable RLOF and CE phases are effectively interchangeable in what follows. In this case, changes in the orbital separation during these phases are typically described by the balance between the orbital energy lost and the binding energy required to eject the envelope (Paczynski Reference Paczynski, Eggleton, Mitton and Whelan1976; Webbink Reference Webbink1984; Iben & Livio Reference Iben1993; Politano Reference Politano2021). This pathway typically produces tighter DWD systems as both episodes of envelope ejection efficiently shrink the orbit. These mechanisms account for the observed populations of close DWDs and their diverse orbital properties (Brown et al. Reference Brown, Gianninas, Kilic, Kenyon and Allende Prieto2016; Noughani et al. Reference Noughani, Nordhaus, Richmond and Wilson2026). The stability of the mass transfer, efficiency of envelope ejection, and initial binary parameters critically determine the final configuration of the system (Olejak et al. Reference Olejak, Belczynski and Ivanova2021). In summary, the DWDs in our simulations arise from a small number of recurring interaction histories, most commonly either two phases of stable RLOF or sequences that include at least one episode of unstable mass transfer leading to a CE (e.g. RLOF+CE, CE+RLOF, or CE+CE). For the purposes of the population analysis below, we therefore group systems into two broad families: a purely stable RLOF channel, comprising binaries that never enter a CE (zero CE episodes), and a CE involving channel, comprising all binaries that experience at least one CE episode, independent of the exact order of stable and unstable mass transfer. Wide DWDs in our models are predominantly produced by the RLOF only channel, whereas the tightest systems almost always originate from CE involving histories.
Another possible pathway for the formation of extremely low-mass (ELM) DWDs involves a two-stage process: the system first undergoes a CE phase triggered by unstable mass transfer, leading to the formation of a compact WD + MS binary. This is then followed by a stable RLOF phase, during which the MS companion transfers mass to the WD, producing an ELM WD (Li et al. Reference Li, Chen, Chen and Han2019). In our simulations, any DWDs formed in this way would belong to the CE involving family, but we do not attempt to isolate or quantify this specific channel separately, any systems formed in this way would simply be included within the broader CE involving population in our analysis.
2.2. Binary population synthesis with COMPAS simulation
To generate synthetic binary populations of WDs, we use compas v02.50.00 (Riley et al. Reference Riley2022; Team COMPAS et al. Reference COMPAS2025), an open-source rapid binary population synthesis (BPS) code. compas simulates the evolution and interaction of millions of stellar binaries using a grid of parametrised physical prescriptions. Similar to other BPS codes such as bse (Hurley et al. Reference Hurley, Tout and Pols2002), STARTRACK (Belczynski et al. Reference Belczynski, Kalogera, Rasio, Taam, Zezas, Bulik, Maccarone and Ivanova2008), BINARY C (Izzard et al. Reference Izzard, Tout, Karakas and Pols2004), SEBA (Toonen, Nelemans, & Portegies Zwart Reference Toonen, Nelemans and Portegies Zwart2012), and COSMIC (Breivik et al. Reference Breivik2020), COMPAS is computationally efficient, enabling rapid exploration of large parameter spaces and statistical studies of binary populations. This efficiency is achieved through simplified evolutionary prescriptions, enabling the evolution of millions of binaries in just a few CPU hours, although at the cost of not resolving the detailed internal structure of individual stars captured by 1D stellar evolution codes such as MESA or other hydrodynamic simulations (Price et al. 2018). While COMPAS is valuable in generating binary populations, it does not resolve detailed stellar structure, hydrodynamics, or dynamical interactions in dense environments, and relativistic effects are included only through approximate orbital decay prescriptions. As a result, predicted binary properties and merger rates are sensitive to model assumptions and should be interpreted statistically rather than as precise predictions for individual systems (Romero-Shaw et al. Reference Romero-Shaw, Hirai, Bahramian, Willcox and Mandel2023; Rauf et al. Reference Rauf, Howlett, Stevenson, Riley and Willcox2024).
While COMPAS was primarily applied to study the formation and evolution of neutron stars and black holes (Broekgaarden et al. Reference Broekgaarden2022; Stevenson & Clarke Reference Stevenson and Clarke2022; Wagg et al. Reference Wagg, Broekgaarden, de Mink, Frankel, van Son and Justham2022), we extend its application to low-mass binaries in this study, focusing specifically on the formation channels and population properties of DWD systems. Our study thus explores the evolution of binaries that produce WDs and establishes a framework for future, more detailed investigations of low and intermediate mass stellar evolution and for calibration against observed DWD populations. Moreover, in the present version of compas, magnetic braking for low-mass stars is not incorporated which mostly acts while low mass stars with convective envelopes are still on the main sequence. This would have made close orbit binaries tighten faster and thereby affecting the formation efficiency and orbital period distribution for low-mass systems (Belloni et al. Reference Belloni, Schreiber, Moe, El-Badry and Shen2024); whereas it has little impact on the wide DWDs.
In this study, we use COMPAS to generate a large population of DWDs by sampling initial binary parameters from observationally motivated distributions. The primary masses are drawn from a Salpeter initial mass function in the range 0.5–15
$\mathrm{M}_{\odot}$
, and the secondary masses are set by a flat initial mass ratio distribution,
$q$
(Salpeter Reference Salpeter1955). This mass range includes primaries that are massive enough to form WDs within a Hubble time while avoiding very low-mass stars that would not evolve off the main sequence. Systems in which either component undergoes core collapse are excluded from our DWD sample. The semi-major axes (a) are sampled from FLATINLOG (i.e. flat in log a) distribution over 0.01–10 AU covering both very close binaries that are almost guaranteed to interact via RLOF and moderately wide systems that can still undergo interaction once one or both stars evolve to giant radii. Much wider binaries, which are unlikely to experience strong interaction or to form close DWDs, are not considered here. Eccentricities are drawn from a thermal distribution of
$f(e)=2e$
with e being the eccentricity (Heggie Reference Heggie1975). All binaries are assumed to have a fixed metallicity (
$Z=0.014$
) throughout the evolution from progenitor binaries to the formation of DWDs. For our main simulation set, we evolve
$10^7$
binaries using the default rapid-synthesis mode of compas, identifying DWDs as systems in which both stars become WDs with He, CO, or ONe-cores (Hurley et al. Reference Hurley, Tout and Pols2002). In compas, WDs are classified as He, CO, or ONe-core objects according to their prior core-burning history (e.g. whether He has been exhausted and whether carbon ignition has occurred), rather than a strict cut on the final WD mass. This is important when interpreting low-mass CO tagged objects and massive CO/ONe remnants in our population. In our simulations, each binary is evolved from zero-age main sequence (ZAMS) up to a maximum age equal to the Hubble time (the default maximum evolution time in compas,
$\sim$
13.8 Gyr), or until both components have become compact remnants. This enables us to analyse their orbital separations, masses, compositions, and mass transfer rates to connect formation channels with observable properties.
Unless otherwise specified, we adopt standard compas physical prescriptions with CE efficiency parameters
$\alpha_{\text{CE}}=1.0$
and
$\lambda_{\text{CE}} = 0.1$
. We do not attempt to recalibrate these parameters to match the observed DWD population. Instead, we treat them as plausible but uncertain choices and also explore how varying
$\lambda_{\text{CE}}$
affects the separation distribution in Section 3. The resulting synthetic population forms the basis of our analysis in subsequent sections.
3. Results
3.1. Effects of evolutionary channels on DWD properties
Our simulations confirm that the evolutionary pathways leading to DWDs can usefully be grouped into two broad categories: (i) systems that evolve solely through stable RLOF and never enter a CE phase, and (ii) systems that experience at least one CE episode. These two interaction channels together account for the large majority of DWDs in our models. In addition to these interaction-driven pathways, our models also produce a small number of DWDs that never undergo sustained RLOF or CE phase. These wide and effectively non-interacting systems are not included in the two formation channels defined above and are not considered hereinafter. Each channel has distinct outcomes on DWD properties, especially in their masses and orbital separation. Figure 1 illustrates this by plotting the final semi-major axis against the individual WD mass formed from each channel. Each point in the plot represents a WD in the binary system, with panel (a) showing the more massive (primary) WD and panel (b) the less massive (secondary) WD in each binary.
Final WD mass as a function of orbital separation for DWDs for
$\alpha_{\text{CE}}=1.0$
and
$\lambda_{\text{CE}} = 0.1$
. Systems evolving solely through stable RLOF are shown in blue, while those involving at least one CE phase are shown in orange. Panel (a) represents the more massive WD and panel (b) represents the less massive companion in each binary system.
Both panels of Figure 1 reveal a bimodal distribution in the orbital separation of DWDs, reflecting their evolutionary histories.
-
• Systems evolving solely through stable RLOF mass transfer tend to remain at wide orbital separations, with final semi-major axes predominantly exceeding
$500\,\mathrm{R}_\odot$
, and in some cases extending beyond
$10^4\,\mathrm{R}_\odot$
, along with a relatively small fraction of RLOF only systems at smaller separations. These close RLOF only binaries arise from net orbital shrinkage during the second stable RLOF episode (see aforementioned discussion of Section 2.1). The component WD masses show a broad distribution, spanning approximately 0.5–1.4
$\mathrm{M}_\odot$
for primary and 0.1–1.15
$\mathrm{M}_\odot$
for secondary component. This reflects the diversity of initial binary parameters and mass transfer histories. A sharp lower mass cut-off for the primary component is evident at approximately
$0.5\,\mathrm{M}_\odot$
, whereas a small gap is present for the secondary component around that value. In the next section, we show that these features arise from the distribution of chemical compositions of WDs. -
• In contrast, systems undergoing at least one CE phase end up in much closer orbits, with their final semi-major axes mainly below
$100\,\mathrm{R}_\odot$
. Their orbital separations are orders of magnitude smaller, while the WD masses are more spread out compared to the stable RLOF only channel. No significant mass gap is observed in this population. These CE involving systems arise from the various RLOF+CE, CE+RLOF, and CE+CE sequences discussed in Section 2, but regardless of the detailed order of events, the presence of at least one CE episode is what drives the strong orbital contraction.
The key feature of both panels in Figure 1 is the well-defined gap in orbital separation between approximately 100 and
$500\,\mathrm{R}_\odot$
. The lack of systems at this intermediate separation in our models highlights a characteristic imprint of binary interaction: the mechanisms that dominate DWD formation preferentially yield either wide or close binaries, with a few systems in between. Stable mass transfer is typically conservative or weakly non-conservative (Soberman, Phinney, & van den Heuvel Reference Soberman, Phinney and van den Heuvel1997), often resulting in orbital expansion, whereas the CE phase is highly efficient at extracting orbital energy, inducing a significant contraction of the orbit. The consistency of this bimodal separation distribution across both primary and secondary WD populations in our simulations strongly suggests that the type of binary interaction (purely stable RLOF versus at least one CE episode) is the dominant factor in setting the final orbital separation of DWD systems. We show later in this section that varying the CE binding parameter
$\lambda_{\text{CE}}$
modifies the depth and width of the gap but not its overall presence.
Current DWD catalogues are majorly incomplete due to strong observational selection effects, particularly for close binaries. However, there is emerging evidence of DWDs at intermediate separations. In particular, Korol et al. (Reference Korol, Hallakoun, Toonen and Karnesis2022)b, Korol, Belokurov, & Toonen (Reference Korol, Belokurov and Toonen2022a) combined population synthesis models with Gaia selected wide DWDs and showed that their best-fitting models require relatively fewer systems at separations near 1 AU in order to match the observed orbital separation distribution. This characteristic scale is consistent with our inference of gap at 100–500
$\mathrm{R}_\odot$
(
$\approx$
0.46–
$2.32$
AU). A detailed, one-to-one comparison with the observed DWD population, including a full analysis of observational selection effects, is beyond the scope of this study. We therefore restrict our discussion to emphasising the order-of-magnitude agreement between our predicted gap and the under-density inferred from Gaia-calibrated Galactic DWD population models relevant for LISA.
The physical origin of the gap is associated with binary interactions during post main-sequence evolution, where the maximum radius (
$R_{\text{max}}$
) of WD progenitors during RGB/AGB evolution can reach upto a few tens to a few hundred solar radii and it eventually sets the characteristic separation scale
$a_{\text{crit}} \simeq R_{\text{max}}/f(q)$
with
$f(q)\equiv R_L/a$
. Systems with
$a\lesssim a_{\text{crit}}$
inevitably undergo RLOF (Lee & Nelson Reference Lee and Nelson1988). If mass transfer is stable and only weakly non-conservative, the orbit widens, yielding final separations
${\gtrsim}500\,\mathrm{R}_\odot$
. A small subset of RLOF-only systems can experience two stable mass transfer phases. During the second episode, mass transfer from the initially less massive star onto the first-formed WD leads to orbital shrinkage, producing relatively close DWDs with separations of a few tens to
$\sim 100 \,\mathrm{R}_\odot$
without a CE phase. In contrast, systems that undergo at least one CE phase experience dynamically unstable mass transfer, resulting in substantial orbital energy loss and final separations
${\lesssim}100\,\mathrm{R}_\odot$
.
3.2. Effects of chemical compositions on DWD properties
To explain the features observed in the orbital separation distributions, especially in tight binaries, we examine the relationship between final orbital parameters and WD core composition. Figure 2 illustrates the final semi-major axis and mass for DWDs categorised by their core composition (He, CO, or ONe). Similar to that in the previous case, panel (a) corresponds to the more massive WD in each binary, while panel (b) shows the less massive component.
Same as Figure 1 except each WD is marked according to their core composition.
Our analysis reveals a strong correlation between core type, mass, and orbital separation.
-
• He-core WDs are confined to low mass ends (
${\lesssim}0.5\,\mathrm{M}_\odot$
), exhibiting a upper mass cut-off. This upper mass limit is consistent with the compas remnant-typing scheme described above, where WDs are classified according to whether their stripped cores ignite He and begin building a CO core. The heavier components in a DWD system are found exclusively in close binaries, originating from a CE evolution. In contrast, secondary He-core WDs appear in both wide and close binaries, indicating diverse formation pathways for the lower-mass component. -
• CO-core WDs dominate the population and span across intermediate to higher mass (
${\gtrsim}0.5\,\mathrm{M}_\odot$
) with a relatively small fraction at lower mass. This is because of compas classification of WDs: once a stripped core ignites He and starts building a CO-core, the remnant is tagged as a CO-core WD even if He burning is incomplete and the total mass remains relatively low. As a result, the region around
$0.5\,\mathrm{M}_\odot$
marks the transition from He-core to CO-core WDs in the models. CO-core WDs are found across the entire spectrum of orbital separations. In wide binaries, the more massive WD is invariably a CO-core WD, as its progenitor must have evolved to at least the tip of the red giant branch (RGB) to ignite He. This behaviour is consistent with standard expectations for He-CO transitions in low and intermediate mass stars, and with previous population synthesis results -
• ONe-core WDs are very few in numbers and occupy the highest mass (
${\gtrsim}1.2\,\mathrm{M}_\odot$
). They are mostly found at intermediate separations and originated from high mass progenitors. In a few systems, the CO-core companion can be comparable in mass or even slightly more massive than the ONe WD, reflecting earlier accretion driven growth of the CO remnant and the fact that the CO-ONe boundary in compas is set by core ignition history rather than a sharp mass threshold.
Figure 2(b) shows the small mass gap for a stable RLOF channel in the secondary WD population (Figure 1b) is associated with the core composition, where He-core WDs exhibit lower masses and CO-core WDs exhibit higher masses. In our models, we do not find any double He-core WD systems in wide binaries formed exclusively through the stable RLOF channel. The formation of such systems requires both stars to undergo envelope stripping during the RGB phase via stable mass transfer, with the evolution precisely stopped prior to He ignition. This process requires a highly specific range of initial stellar masses, orbital separations, and mass transfer efficiencies. In practice, the secondary star typically evolves to a more advanced evolutionary stage, resulting in a CO-core WD, thereby preventing the formation of double He-core wide binaries. The absence of such systems in our population should therefore be interpreted as a consequence of our adopted stability criteria and mass transfer prescriptions, rather than as a strict prediction that wide double He-core DWDs cannot form in nature.
3.3. Effects of mass transfer rates on DWD properties
To further understand the features observed in the orbital separation distributions, especially in tight binaries, we investigate the role of mass transfer efficiency. In particular, we examine how the efficiency of stable mass transfer and the associated angular-momentum loss map onto the final orbital separations, which is expected to produce a correlation between mass transfer rate and post-interaction separation in standard binary-evolution theory. Figure 3 shows the final semi-major axis and mass for DWDs, with the colour bar indicating the mass transfer rate during the most recent stable RLOF episode. Note that for this analysis, we restrict the sample to systems that underwent at least one stable RLOF mass transfer phase. This removes many WDs including the primary components with masses
${\lesssim}0.5\,\mathrm{M}_\odot$
, which are typically formed without any stable RLOF episode.
Same as Figure 1 except that the donor mass transfer rates during the last stable mass transfer episode is shown in the colour bar.
The key observations from the plots are as follows.
-
• Wide binaries are characterised by the highest mass transfer rates. This is a direct consequence of the compas prescription for stable mass transfer. When a star overflows its Roche lobe, compas initiates stable mass transfer on a thermal (or, for some systems, nuclear) timescale. In these wide systems, the mass transfer is often nearly conservative, brief, and intense, resulting in high donor mass transfer rates. This interaction does not significantly alter the orbit, allowing both stars to subsequently evolve independently and form a wide DWD.
-
• Close and interacting binaries, on the other hand, can evolve into DWDs by mass transfers with lower to higher mass transfer rates. Comparing Figure 3 with Figure 2, it is evident that for close DWDs with the primary component being a CO-core WD and the secondary being a He-core WD, the mass transfer rate is the lowest. For double CO-core WD in close binaries, the mass transfer rate is relatively higher, and it increases as the final orbital separation decreases. This anti-correlation between final separation and mass transfer rate is consistent with the expectation from angular-momentum conservation: phases of efficient angular-momentum loss, whether through non-conservative mass loss or the onset of a CE, drive stronger orbital decay and produce tighter systems. In our models, close double CO-core WDs tend to descend from such high mass transfer, high angular-momentum loss episodes. No strong correlation is found between WD mass and mass transfer rate.
Together, the interconnection between mass transfer rate and orbital separation provides a quantitative illustration of how our adopted mass transfer and angular-momentum loss prescriptions shape the final properties of DWDs, and is consistent with the general expectations from angular-momentum conservation in interacting binaries.
3.4. Dependence on the CE binding parameter
The characteristics of the DWD population created by our models are determined primarily by the selected CE parameters, specifically by the product
$\alpha_{\text{CE}}\lambda_{\text{CE}}$
, which controls how much orbital energy is available to unbind the envelope. In our fiducial model, we use
$\alpha_{\text{CE}}=1$
and
$\lambda_{\text{CE}}=0.1$
, corresponding to relatively tightly bound envelopes. To assess the robustness of our main outcomes, we perform an additional simulation where we keep
$\alpha_{\text{CE}}=1$
but increase the envelope-structure parameter
$\lambda_{\text{CE}}=0.5$
. This choice effectively increases the product
$\alpha_\mathrm{CE}\lambda_\mathrm{CE}$
, decreases the magnitude of the binding energy of the envelope, allowing a larger fraction of binaries to survive CE ejection and form DWDs. In contrast, the
$\lambda_{\text{CE}}$
parameter does not affect the stability of RLOF mass transfer which instead depends on the donor’s evolutionary stage along with mass ratio, the orbital separation, Roche lobe size, and the adopted mass and angular-momentum prescriptions. Consequently, this model produces a significantly larger number of CE-involving DWDs, while the number of DWDs formed exclusively through the RLOF-only channel remains unchanged relative to the aforementioned fiducial model.
Figure 4 illustrates the orbital separations obtained at the end of the RLOF-only and CE-involving channels in the two CE models. The distribution of separations remains clearly bimodal in both cases, with systems that have experienced only RLOF occupying mostly the wide orbits and the CE-involving systems making up the close DWDs. The lack of systems at intermediate separations between tens and a few hundred
$\text{R}_\odot$
is intact, although its depth and exact boundaries shift when
$\lambda_{\text{CE}}$
is altered. Increasing
$\lambda_{\text{CE}}$
allows CE survivors to emerge at somewhat larger separations, thus partially filling the gap, which is anticipated from the fact that the orbital shrinkage is lesser when more energy is being used to unbind the envelope.
Same as Figure 1 except for
$\alpha_{\text{CE}}=1$
and
$\lambda_{\text{CE}}=0.5$
.
Overall, the presence of a separation gap in DWD population is a robust feature of our models, while the precise appearance of the gap (its width, depth, and the relative numbers of systems on either side) does depend on the underlying CE parameters. In the remainder of the paper, we therefore interpret the separation gap as a qualitative prediction of our binary-evolution framework, and regard the precise shape of the distribution as a useful diagnostic for future calibration against observed DWD populations.
3.5. Illustrative evolutionary histories
In order to verify the trends seen in the mass vs. semi-major axis distributions are genuine physical effects and not numerical artifacts, we inspect full evolutionary histories of selected systems drawn from our simulated population. Below, we present two cases that detail the previously mentioned channels.
The first case involves a pair of close DWD, which is composed of a low mass CO WD with a more massive He-core companion. The system originates from an initially unequal mass main sequence binary. The initially more massive star is the first to fill its Roche lobe as a subgiant and undergoes stable mass transfer, being stripped down to a low mass naked He star while while the companion that accretes mass to become the more massive component. At this stage, the orbit widens and circularises. When the companion (now more massive) later ascends the giant branch, its RLOF is dynamically unstable and causes CE to occur. The envelope is ejected, leaving behind a close binary comprising the He star and a newly formed He-core WD. The He star then completes core He burning and cools as a WD. As its stripped core has ignited He and begun to build a CO-core, compas classifies the remnant as a CO-core WD even though its mass is still low. This example demonstrates how early envelope stripping and incomplete He burning naturally yield low mass CO tagged (hybrid He-CO) WDs in close binaries that have come about through a combination of stable RLOF and CE evolution.
The second case is a close and massive DWD composed of a high-mass CO-core WD and a slightly less massive ONe-core star. In this scenario, the star which is initially the less massive becomes a strong mass gainer during an early phase of stable RLOF, growing into a massive stripped He star. After He exhaustion, this star exposes a massive CO-core and cools as a CO WD. The evolution of the companion into a massive giant occurs much later, which results in the initiation of unstable mass transfer leading to a CE phase. The envelope is ejected, and the core of the giant, which has undergone carbon ignition, is left behind as an ONe-core WD in a close orbit around the pre-existing massive CO WD. In this scenario, the CO-core WD is just a bit heavier than its ONe-core partner due to the prior accretion driven growth and the fact that the CO-ONe boundary in compas is determined by core ignition history rather than a sharp mass threshold.
These illustrative full evolutionary histories confirm that the apparently unusual regions of parameter space in our population, such as low mass CO WDs and CO-ONe pairs where the CO component is the more massive star, arise as natural outcomes of well-defined mass transfer and CE channels in the binary evolution model. Furthermore, the core type classification method that is used in COMPAS must be taken into account in when analysing the mass and composition distributions of synthetic DWD populations as it is based on the prior burning history of the stripped core.
4. Discussion
In this study, we use the compas simulation to perform a comprehensive population synthesis of how DWDs form and evolve. Our population synthesis results demonstrate that the evolutionary pathway of a binary system shapes the final properties of DWD pairs. A similar low density in the separation distribution of DWDs has been reported previously in combined population synthesis and Gaia analyses, and our models reproduce this feature and quantify its dependence on binary interaction physics. The most salient feature in our simulations is the significant bimodal distribution in orbital separation, characterised by a clear gap in 100–500
$\mathrm{R}_\odot$
. This feature emerges directly from the underlying physics of mass transfer: a stable RLOF typically preserves or expands the orbit, whereas CE evolution induces substantial orbital contraction. The way this gap splits the close and wide DWD populations provides a powerful, testable signature for future observatories. GW detectors like LISA will be sensitive to the close binary population, while high-precision astrometric surveys (e.g. Gaia, Roman) can characterise wide binaries, enabling stringent tests of binary evolution models (Nelemans et al. Reference Nelemans, Portegies Zwart, Verbunt and Yungelson2001b; Toonen et al. Reference Toonen, Nelemans and Portegies Zwart2012). In particular, we find that varying the CE binding parameter
$\lambda_{\text{CE}}$
changes the depth and width of the separation gap but not its overall presence, indicating that the gap is a robust outcome of our models while its morphology encodes details of the CE physics.
The core composition of WDs further reflects their evolutionary pathways. Double He-core WDs, with masses
${\lesssim}0.5\,\mathrm{M}_\odot$
, in our models form only in close binaries, as they require early envelope stripping via at least one CE phase that also shrinks the orbit. In contrast, CO- and ONe-core WDs are massive and form across a much wider range of orbital separations, reflecting a more diverse set of evolutionary channels (Iben & Tutukov Reference Iben1985; Han Reference Han1998; Toonen et al. Reference Toonen, Hollands, Ga¨nsicke and Boekholt2017). Notably, our simulation does not produce wide double He-core WD systems. The formation of such a pair would require two consecutive episodes of stable mass transfer that precisely truncate the evolution on RGB for both stars, a scenario demanding a finely tuned set of initial conditions that is statistically improbable. In our simulations such finely tuned configurations do not occur, and we therefore do not find any wide double He-core DWDs. This absence should, however, be interpreted as a consequence of our adopted stability criteria and mass transfer prescriptions rather than as proof that such systems cannot exist in nature.
The observed correlation between the mass transfer rate and final orbital separation offers critical insight into mass transfer physics. High mass transfer rates are associated with brief conservative mass transfer events that result in wide orbits. Conversely, within the population of close binaries, we find that orbital separation decreases with increasing mass transfer rates. Within the population of close binaries we find that systems with smaller final separations tend to have experienced more intense stable mass transfer episodes shortly before DWD formation. This trend is consistent with the expectation that phases of efficient angular-momentum loss during mass transfer, whether through non-conservative mass loss, enhanced winds, or the onset of a CE drive stronger orbital decay and produce tighter DWDs.
One particularly interesting outcome is the substantial population of wide DWDs produced in our simulations. Although their mutual gravitational interaction is weak, these systems remain bound and in principle can serve as natural laboratories for testing fundamental physics. Recent studies incorporating Gaia data have shown that wide binaries are excellent probes for modified gravity theories, such as MOND, providing complementary constraints to galaxy dynamics (Pittordis & Sutherland Reference Pittordis and Sutherland2019; Banik et al. Reference Banik, Pittordis, Sutherland, Famaey, Ibata, Mieske and Zhao2024).
Several important assumptions directly affect the properties of the simulated DWD population. The evolution of the CE phase plays a major role: in the standard energy formalism the product
$\alpha_{\text{CE}}\lambda_{\text{CE}}$
controls how much orbital energy is available to unbind the envelope. Higher values of
$\alpha_{\text{CE}}\lambda_{\text{CE}}$
make envelope ejection easier, increasing the overall number of CE survivors and favouring somewhat wider post-CE systems, whereas lower values lead to stronger orbital contraction and a higher rate of stellar mergers but also enhance the formation of the tightest DWDs that contribute to the LISA band. The stability criterion for RLOF is equally important. A more stringent criterion causes the systems to undergo CE evolution, leading to the formation of tight binaries. In contrast, a more relaxed criterion allows stable RLOF and increases the number of wide binaries. Since the stability boundary is still poorly constrained, variations in this criterion can significantly shift the relative contributions of the RLOF only and CE involving channels in any given population synthesis model. The efficiency of mass transfer (
$\beta$
, the fraction of transferred mass that is accreted by the companion) also influences the outcomes. In many commonly used prescriptions, nearly conservative transfer (
$\beta = 1$
) tends to favour orbital expansion for mass ratios close to unity, whereas non-conservative transfer (
$\beta\lt1$
) can remove substantial orbital angular momentum through ejected material and often leads to orbital contraction, especially for unequal mass systems.
In addition to these binary interaction parameters, outputs are also affected by global parameters including the metallicity and the star formation rate (SFR). Lower metallicity environments typically form more massive WDs and alter CE survival rates, whereas higher metallicity favours low mass remnants and wider systems. Similarly, a constant SFR assumption provides a time averaged view of the Galactic population, but a bursty or declining SFR would change the relative proportions of young, close binaries to older, wider binaries. For simplicity, we fixed the metallicity and SFR in this work, and a detailed systematic study of their effects is reserved for future work. Nonetheless, within the range of parameter variations we have explored, the qualitative features of the population, such as the presence of a separation gap between RLOF only and CE involving systems and the broad mass and core type trends remain robust, even though the detailed morphology and relative weights of the different channels can change significantly.
4.1. DWD formation efficiency and comparison with previous studies
In our fiducial COMPAS run (
$\alpha_{\text{CE}}\lambda_{\text{CE}}=0.1$
), we evolve
$N_\mathrm{bin} = 10^{7}$
ZAMS binaries with
$Z=0.014$
and a total initial stellar mass
$M_\mathrm{init,tot} = 2.021045\times 10^{7}\,\mathrm{M}_\odot$
. Among them,
$N_\mathrm{ DWD} = 329\,147$
systems form DWDs, corresponding to a formation efficiency
\begin{equation}f_\mathrm{DWD} = \frac{N_\mathrm{DWD}}{N_\mathrm{bin}} \simeq 3.29\times 10^{-2},\end{equation}
and a DWD yield per unit stellar mass formed
\begin{equation}\eta_\mathrm{DWD} = \frac{N_\mathrm{DWD}}{M_\mathrm{init,tot}} \simeq 1.63\times 10^{-2}\,\mathrm{M}_\odot^{-1}.\end{equation}
Breaking this down by interaction channel,
$N_\mathrm{DWD,RLOF} = 176\,915$
systems (
$\simeq 54\%$
) form through only stable RLOF, while
$N_\mathrm{DWD,CE} = 115\,182$
systems (
$\simeq 35\%$
) experience at least one CE episode. The residual
$N_\mathrm{DWD,other} = 37\,050$
systems (
$\simeq 11\%$
) form without any mass transfer episode.
Assuming the birth rate of DWDs (
$\mathcal{R}_\mathrm{DWD}$
) scales with the SFR of Milky Way (MW), it can be expressed as
\begin{equation}\mathcal{R}_\mathrm{DWD} = \eta_\mathrm{DWD}\times\text{SFR}_{\text{MW}}.\end{equation}
Following the recent estimate by Siegert et al. (Reference Siegert, Pleintinger, Diehl, Krause, Greiner and Weinberger2023), we choose
$\text{SFR}_{\text{MW}} = 6\pm2 \,\mathrm{M}_\odot\,{\rm yr}^{-1}$
, which yields
$\mathcal{R}_\mathrm{DWD} \simeq (9.6\pm3.2)\times$
$10^{-2}\,\mathrm{yr}^{-1}$
. With the channel splits, these values turn out to be
$\mathcal{R}_\mathrm{DWD,RLOF}\simeq (5.28\pm1.76)\times 10^{-2} \,\mathrm{yr}^{-1}$
,
$\mathcal{R}_\mathrm{DWD,CE} \simeq (3.42\pm1.14)\times 10^{-2} \,\mathrm{yr}^{-1}$
, and
$\mathcal{R}_\mathrm{DWD,other} \simeq (1.08\pm0.36)\times 10^{-2} \,\mathrm{yr}^{-1}$
. In the second simulation with
$\alpha_{\text{CE}}\lambda_{\text{CE}}=0.5$
, these numbers change to
$\mathcal{R}_\mathrm{DWD,RLOF}\simeq (5.28\pm1.76)\times 10^{-2} \,\mathrm{yr}^{-1}$
,
$\mathcal{R}_\mathrm{DWD,CE} \simeq$
$ (9.47\pm3.16)\times 10^{-2} \,\mathrm{yr}^{-1}$
, and
$\mathcal{R}_\mathrm{DWD,other} \simeq (1.08\pm0.36)\times 10^{-2} \,\mathrm{ yr}^{-1}$
. This simple scaling can be compared with earlier Galactic DWD population synthesis studies. Using SeBa binary population synthesis code, Nelemans et al. (Reference Nelemans, Yungelson, Portegies Zwart and Verbunt2001a) reported a Galactic birth rate of close DWDs of
$\simeq 5\times 10^{-2}\,\mathrm{ yr}^{-1}$
with a total present day population of
$\simeq 2.5\times 10^{8}$
systems. Other rapid binary-evolution models that adopt the prescriptions of Hurley et al. (Reference Hurley, Tout and Pols2002), typically predict current DWD birth rates of
$3.21\times 10^{-2}\,\mathrm{yr}^{-1}$
in the thin disc, depending on the assumed star-formation history (Yu & Jeffery Reference Yu and Jeffery2010). Despite subtle differences in input physics and calibration, our compas-based rates are thus broadly consistent with previous population synthesis results obtained using independent codes and assumptions within the expected factor of few systematic uncertainties.
5. Conclusions
In this work, we have used the rapid binary population synthesis code COMPAS to explore how evolutionary pathways, WD core compositions, and mass transfer rates jointly shape the final properties of DWD systems. Our models reproduce the bimodal separation distribution and intermediate under-density previously identified in studies combining population synthesis with Gaia data. Thereby we show that this feature arises naturally from the coexistence of two broad families of interaction channels: wide systems that evolve purely through stable RLOF and close systems that experience at least one CE phase. Within this unified framework we quantify how the positions of DWDs in the mass–separation plane correlate with their core types (He, CO, ONe) and with the intensity of their last stable mass transfer episode, and we demonstrate that the presence of a separation gap between wide RLOF only and close CE involving systems is robust to reasonable changes in the CE structure parameter, even though its detailed morphology does depend on the adopted CE physics.
These results provide a physically interpretable map between the observed parameters of DWDs (masses, separations, and core types) and the underlying interaction histories predicted by COMPAS like binary evolution modes. They therefore offer a concrete framework for using future observations to calibrate uncertain aspects of binary physics, such as the stability of mass transfer and the efficiency of CE ejection. Looking ahead, forthcoming observatories will greatly improve our empirical picture of the DWD population: space-based GW missions such as LISA will systematically survey the Galactic population of close, interacting binaries, while astrometric missions such as Gaia and Roman will enhance the detection and characterisation of wide DWDs. Combining these data sets with population synthesis models like those presented here will be crucial for refining binary evolution prescriptions and may eventually allow tests of gravitational physics on binary scales.
Acknowledgements
The authors would like to thank the anonymous reviewer for their insightful comments, which have helped to improve the manuscript content. We thank T. Bulik and D. Rosińska from University of Warsaw for useful discussion during the compilation of the work. SK acknowledge funding from the National Center for Science, Poland (grant no. 2023/49/B/ST9/02777).
Data availability statement
The data used in this article are obtained from the compas v02.50.00 simulation, which will be shared on reasonable request to the corresponding author.
Open access
