© The Authors 2024

1 Introduction

The Milky Way contains ~10⁷ double white dwarf (DWD) binaries, which emit nearly monochromatic gravitational waves (GWs) within the mHz band (Breivik et al. 2018; Nelemans et al. 2001b; Ruiter et al. 2010; Lamberts et al. 2019). The Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017), which is scheduled for launch in 2035, will complement electromagnetic instruments with gravitational wave (GW) detection, and will thus provide a unique opportunity to study this largely unobserved population. Unresolvable DWDs are expected to produce stochastic noise that will exceed LISA’s instrumental noise in the ~0.1–2 mHz range. Additionally, up to tens of thousands of binaries are expected to be individually resolved. This makes the analysis of LISA data challenging. This challenge has spurred the development of Global Fit pipelines (Littenberg & Cornish 2023), which are designed to infer this persistent yet unknown signal present in the data. The detection and characterisation of this population will help us understand the formation and evolution of white dwarfs (WDs), which is one of the eight objectives of the LISA mission (Colpi et al. 2024). Now that LISA has been adopted by the European Space Agency, it is timely to construct the tools for the analysis of LISA data and its astrophysical interpretation. Similarly to the recent study of Scaringi et al. (2023), which highlighted observational signatures of cataclysmic variables within LISA’s confusion noise, our work anticipates distinct observational features, attributable in our case to tidal effects in mass-transferring binaries.

The evolution of DWDs is characterised by various physical processes that can significantly influence their behaviour and observational signatures. One of the key aspects is the effect of matter interactions, which becomes relevant early in the evolution, due to the reduced compactness of WDs compared to neutron stars and black holes. Similar to binary stars, tidal torques work to synchronise the rotation period of the WDs with the orbital period of the binary (Zahn 1977; Verbunt & Hut 1983). While expressions for the intensity of these torques within DWDs have been proposed (Campbell 1984), they are generally valid only under the assumption of low asynchronicity, and the precise intensity of these torques remains uncertain. As the separation between the WDs decreases, the lighter WD may overflow its Roche lobe, leading to mass transfer onto the heavier WD. This scenario can result in one of two outcomes. If the mass transfer is dynamically unstable, the binary may merge rapidly, potentially triggering a Type Ia supernova (Webbink 1984; Woosley & Weaver 1986; Fink et al. 2007, 2010; Guillochon et al. 2010; Kromer et al. 2010; Maoz et al. 2014; Shen et al. 2018). If instead the mass transfer is stable, then an AM CVn system can be formed (Nather et al. 1981; Nelemans et al. 2001c; Solheim 2010; Tutukov & Yungelson 2002). These stable interacting binaries offer unique opportunities for multi-messenger observations (Nelemans et al. 2004; Shah & Nelemans 2014; Kupfer et al. 2019, 2023).

Marsh et al. (2004) proposed the first semi-analytical model for the evolution of DWDs including tidal torques and mass transfer, in addition to GWs. This approach was extended in Gokhale et al. (2007) to include the evolution of the lighter WD spin. Finally, Kremer et al. (2015) proposed a variation of the previous models with a more refined treatment of mass transfer based on the ballistic accretion of particles (Sepinsky & Kalogera 2014). Those works highlighted the effect of tidal torques on the evolution of DWDs, and in particular that strong tides can help stabilise the dynamics of mass transfer, preventing rapid mergers. We note that in these analyses, the stability of mass transfer is determined using a simple and thus potentially unrealistic criterion based on the rate of mass lost by the donor.

Systems surviving mass transfer are predicted to outspiral (i.e. they will have a negative orbital frequency derivative). This distinctive signature can be used to identify the subpopulation of mass-transferring DWDs with LISA, as was proposed in Breivik et al. (2018), Kremer et al. (2017), Biscoveanu et al. (2023) and Yi et al. (2024). In particular, in Biscoveanu et al. (2023), the authors simulate a population of DWDs until the onset of mass transfer using the COSMIC population synthesis code (Breivik et al. 2020) and use the semi-analytical model to compute the long-term evolution of DWDs. They explore three distinct assumptions regarding the intensity of tidal coupling, and perform a population analysis on these DWDs to demonstrate the potential of LISA observations to probe the intensity of tidal coupling within these systems. They also investigate the impact of different assumptions on the common envelope phase during the evolution of DWDs progenitors, complementing the work of Kremer et al. (2017). Those works shed light on the potential of LISA to inform us about the outcomes of stellar evolution and about the properties of DWDs. However, the tools currently used to perform population analyses (Mandel et al. 2019) are not applicable to DWDs because the signals overlap and, therefore, are not statistically independent. Moreover, it is necessary to account for the dependence of the contribution of DWDs to the confusion noise on the properties of the population that we seek to infer.

Alternatively, approaches to measure mass transfer effects in individual binaries have been developed (Breivik et al. 2018; Yi et al. 2024). For non-accreting systems, it has been suggested that it would be possible to measure the chirp mass, and distinguish between WDs and black holes in Galactic binaries (Sesana et al. 2020; Sberna et al. 2021), to measure tidal effects (Shah et al. 2015), and possibly even to perform tests of general relativity (GR) (Littenberg & Yunes 2019).

In our study, we investigate for the first time the impact of tidal effects within DWDs on the confusion noise that they produce and on the properties of resolvable sources. Starting from a population of DWDs at formation predicted by the population synthesis code of Lamberts et al. (2019), we use a semi-analytical model, detailed in Sect. 2, similar to those of Marsh et al. (2004) and Sberna et al. (2021) to evolve the population until today. We discuss the main features of the long-term evolution of accreting systems in Sect. 3. In this semi-analytical model, tidal torques depend on (i) the spin of the WDs and (ii) the tidal synchronization timescale. The latter measures the efficiency of tidal coupling: the shorter the synchronization timescale, the stronger the tidal interaction. We introduce a universal parametrisation for the tidal synchronization timescale, which allows the intensity of tidal torques to depend on a single parameter for all DWDs (apart from the spin). We generate mock populations exploring a total of 30 combinations of initial spin and tidal synchronization timescale in Sect. 4, and in Sect. 5, we estimate the confusion noise due to DWDs along with the number of detectable binaries, using the methodology of Karnesis et al. (2021), Timpano et al. (2006), Crowder & Cornish (2007), Nissanke et al. (2012) and Lamberts et al. (2019).

We find significant impacts on the DWD population. In particular, if WDs are born slowly spinning and if tidal effects are strong, we anticipate a detectable excess of events around 1 mHz due to accreting binaries that outspiral and accumulate around that frequency. This excess affects both the confusion noise and the resolvable binaries. The latter is in agreement with the findings of Biscoveanu et al. (2023). Moreover, we estimate how many resolvable systems have measurable GW-frequency derivative. We find astrophysical effects in the evolution (i.e. other than GW radiation) to be relevant for all systems with measurable GW-frequency derivative, even for non-accreting systems. Thus, measuring properties of the binary without accounting for effects other than GW radiation will lead to biased estimates. Finally, in Sect. 6, we assess the robustness of our conclusions against systematic effects by considering modifications to our semi-analytical model. We find our results to be generally robust, with the main uncertainty lying in the validity of the criterion we employ to determine the stability of mass-transferring DWDs. We present our concluding remarks in Sect. 7 and give some extra details on oscillatory behaviours in the evolution of mass-transferring DWDs in Appendix A.

2 Evolution of double white dwarf binaries

Our DWD population is based on the catalogue created by Lamberts et al. (2019). This model uses a zoom-in simulation of a galaxy with properties similar to the Milky Way to provide a metallicity-dependent star formation rate that naturally reproduces a thin and thick galactic disk, as well as a bulge and a halo (Wetzel et al. 2016). This was combined with the output of the binary population synthesis code BSE (Hurley et al. 2002) with standard parameters. The population was evolved until the formation of the DWD, and we only kept the systems that have reached that stage by the present day. From this point on, we used the semi-analytical described below to evolve the DWDs until today. In comparison, Kremer et al. (2017) and Biscoveanu et al. (2023) evolved DWDs until the onset of mass transfer with a population synthesis code and, from there, used a semi-analytical model similar to ours to evolve the binaries for 10 Gyr. Our model naturally produces combinations of WDs with different core compositions (helium, carbon/oxygen or neon), which have different evolutionary timescales and typical masses and sky localisations. In the remainder of this paper, we consider the full population of systems, but do not distinguish the different subtypes of systems. In Sect. 6, we discuss the impact of the assumptions we make to model the evolution of DWDs. We use m_i, ω_i, J_i and R_i (i = 1, 2) to represent respectively the masses, spins, angular momenta and radii of the WDs, with m₁ ≥ m₂. We define the mass ratio q = m₂/m₁ ≤ 1. We denote by M_t = m₁ + m₂ the total mass of the binary, a its separation, ω_o its orbital frequency and J_o its orbital angular momentum. We assume all binaries to be quasi-circular.

We begin by specifying our model for mass transfer. Using R_L for the Roche lobe radius of the lighter WD, we define the Roche lobe overfill factor as $$\mathrm{\Delta }=R_2-R_L.$$(1)

We assume that mass is exchanged when the lighter WD overflows its Roche lobe (i.e. when Δ > 0). We use Eggleton’s approximation (Eggleton 1983) for R_L, and following Marsh et al. (2004), we take $${\dot{m}}_2=\{\begin{array}{l}-{\mathcal{F}}_{\text{acc }}(m_2,m_1,a,R_2){\mathrm{\Delta }}^3,\text{ if }\mathrm{\Delta }>0,\\ 0\text{ otherwise, }\end{array}$$(2)

where ℱ_acc is defined in Eq. (10) of Marsh et al. (2004). We limit the rate of mass accreted by the heavier binary to the Eddington limit or below, using the approximation of Han & Webbink (1999) for the latter: $${\dot{m}}_1=min(-{\dot{\mathrm{m}}}_2,{\dot{\mathrm{m}}}_{\mathrm{E}\mathrm{d}\mathrm{d},1}).$$(3)

Therefore, part of the mass might be lost during mass transfer episodes, i.e. we do not necessarily have ${\dot{M}}_t=0$.

The evolution of the binary is governed by the angular momentum balance equation. We assume that the system loses angular momentum only via GWs and loss of matter: $${\dot{J}}_o+{\dot{J}}_1+{\dot{J}}_2=-{\dot{J}}_{\mathrm{G}\mathrm{W}}-{\dot{J}}_{\text{loss}}.$$(4)

The angular momentum of a WD is related to its spin through $$J_i=k_i{m}_i{R}_i{\omega }_i,$$(5)

where k_i is a numerical factor such that $k_i{m}_i{R}_i^2$ is the moment of inertia of the WD. We use the fit from Marsh et al. (2004): k_i = 0.1939(1.44885 − m_i)^0.1917. We assume that the angular momenta of the WDs evolve due to mass transfer and tidal torques as $${\dot{J}}_i=j_i{\dot{m}}_i-\frac{k_i{m}_i{R}_i^2}{{\tau }_{s,i}}({\omega }_i-{\omega }_o).$$(6)

In the above equation τ_s,i is the tidal synchronization timescale. The smaller the timescale, the more efficient tidal torques are in synchronising the spin of the WD with the orbital frequency. In Yu et al. (2020), the authors highlighted that the linear approximation for tidal torques in Eq. (6) is not expected to hold for the higher frequency DWDs observable by LISA. However, in this work, we adhere to this linear approximation to conduct an initial assessment of the impact of tidal effects on the population of DWDs. While the heavier WD might be spun up through accretion from the companion star that will form the lighter WD, the latter is expected to be slowly rotating at birth, and so its specific angular momentum should contribute little to the evolution of the DWD. Thus, for simplicity, we take j₂ = 0. However, for completeness, as described below, we consider the impact of the donor being rapidly rotating at birth, and, in Sect. 6, we explore another possibility for j₂. For the angular momentum of the more massive WD we use the model of Verbunt & Rappaport (1988): $$j_1=\sqrt{G{m}_1{R}_h}.$$(7)

Here, R_h is the radius of the orbit of the matter around the accreting WD. A fit for R_h is provided in Verbunt & Rappaport (1988). This formula corresponds to direct impact accretion. After some time, an accretion disk might form. This occurs when the minimum radius reached by the stream, R_min, exceeds the radius of the accretor. In this case, R_h is replaced by R₁ in Eq. (7). A fit to the minimum radius of the stream is given by Nelemans et al. (2001a): $$\begin{array}{rl}\frac{R_{min}}{a}=& 0.04948-0.03815\text{\hspace{0.17em}log}(q)+0.04752{\text{\hspace{0.17em}log}}^2(q)\\ & -0.006973{\text{\hspace{0.17em}log}}^3(q).\end{array}$$(8)

The evolution equation for ω_i is obtained by taking the derivative of Eq. (5), $${\dot{J}}_i=k_i{m}_i{R}_i^2{\omega }_i({\lambda }_i\frac{{\dot{m}}_i}{m_i}+\frac{{\dot{\omega }}_i}{{\omega }_i}),$$(9)

with ${\lambda }_i=1+2\frac{\mathrm{d}\text{\hspace{0.17em}log}\left(R_i\right)}{\mathrm{d}\text{\hspace{0.17em}log}\left(m_i\right)}+\frac{\mathrm{d}\text{\hspace{0.17em}log}\left(k_i\right)}{\mathrm{d}\text{\hspace{0.17em}log}\left(m_i\right)}$, and equating it with Eq. (6), which yields $${\dot{\omega }}_i=(\frac{j_i}{k_i{R}_i^2}-{\lambda }_i{\omega }_i)\frac{{\dot{m}}_i}{m_i}-\frac{{\omega }_i-{\omega }_o}{{\tau }_s}.$$(10)

Following Marsh et al. (2004), we limit the WD spin to the breakup rate of the WD defined as $${\omega }_{i,K}=\sqrt{\frac{G{m}_i}{R_i^3}}.$$(11)

If the WD spin reaches this value, then we enforce ω_i = ω_i,K and ${\dot{\omega }}_i={\dot{\omega }}_{i,K}$. The latter is achieved by changing the tidal synchronization timescale so that Eq. (10) yields the desired evolution equation for ω_i.

To obtain the evolution equation for the separation, we start from Kepler’s law ${\omega }_o=\frac{\sqrt{Gm}}{a^{3/2}}$, so that the orbital angular momentum can then be written $J_o=\sqrt{\frac{Ga}{M_t}}{m}_1{m}_2$. Therefore, $$\frac{{\dot{J}}_o}{J_o}=\frac{1}{2}\frac{\dot{a}}{a}+\frac{{\dot{m}}_1}{m_1}+\frac{{\dot{m}}_2}{m_2}-\frac{1}{2}\frac{{\dot{M}}_t}{M_t}.$$(12)

At Newtonian order: $$j_{\mathrm{G}\mathrm{W}}=\frac{32}{5}\frac{G^3}{c^5}\frac{m_1{m}_{2}m}{a^4}{J}_o.$$(13)

Following van Haaften et al. (2012), we assume an isotropic reemission of the lost mass and take $${\dot{J}}_{\text{loss }}=q\frac{{\dot{M}}_t}{M_t}{J}_o.$$(14)

Finally, putting together Eqs. (4), (6), (12), (13), and (14), we get $$\begin{array}{rl}\frac{\dot{a}}{2a}=& -\frac{32}{5}\frac{G^3}{c^5}\frac{m_1{m}_2{M}_t}{a^4}-\frac{{\dot{m}}_1}{m_1}(\frac{m_1{j}_1}{J_{\text{orb }}}+1)-\frac{{\dot{m}}_2}{m_2}\\ & +\frac{k_1{m}_1{R}_1^2}{J_{\text{orb }}{\tau }_{s,1}}({\omega }_1-{\omega }_o)+\frac{k_2{m}_2{R}_2^2}{J_{\text{orb }}{\tau }_{s,2}}({\omega }_2-{\omega }_o)+(\frac{1}{2}-q)\frac{{\dot{M}}_t}{M_t}.\end{array}$$(15)

We obtain the corresponding equation for the GW frequency at which the binary emits, combining $f_{\mathrm{G}\mathrm{W}}=\frac{{\omega }_o}{\pi }$ with Eq. (15): $$\begin{array}{rl}\frac{{\dot{f}}_{\mathrm{G}\mathrm{W}}}{f_{\mathrm{G}\mathrm{W}}}& =\frac{96}{5}\frac{G^3}{c^5}\frac{m_1{m}_{2}m}{a^4}+\frac{3{\dot{m}}_1}{m_1}(\frac{m_1{j}_1}{J_o}+1)+\frac{3{\dot{m}}_2}{m_2}\\ & -\frac{3k_1{m}_1{R}_1^2}{J_o{\tau }_{s,1}}({\omega }_1-{\omega }_o)-\frac{3k_2{m}_2{R}_2^2}{J_{\mathrm{o}\mathrm{r}\mathrm{b}}{\tau }_{s,2}}({\omega }_2-{\omega }_o)+(3q-1)\frac{{\dot{M}}_t}{M_t}.\end{array}$$(16)

Finally, we need to specify a mass-radius relation for WDs. We take the commonly used model by Eggleton, defined in Eq. (15) of Verbunt & Rappaport (1988). We note that for small mass WDs (≲0.01 M_⊙) finite-entropy effects are expected to become relevant (Deloye & Bildsten 2003). The systems we evolve do not reach this regime, so we neglect such effects. Thus defined, our model has two sets of unspecified parameters: the tidal synchronization timescale of WDs (τ_s,i) and the spin of WDs after formation (ω_i,0). Indeed, our population synthesis code predicts the white dwarf masses and separation at the formation of the binary, but not their spins.

As a general rule, Campbell (1984) found that the tidal synchronization timescale in a binary scales as: $${\tau }_{s,i}\propto {\left(\frac{m_i}{m_{-i}}\right)}^2{\left(\frac{a}{R_i}\right)}^6,$$(17)

where m_−i is the mass of the other binary component. Thus, we define a reference value τ_s,ref for the heavier WD in a reference system with m₁ = 0.6 M_⊙ and m₂ = 0.25 M_⊙ at the moment it overfills its Roche lobe, and use the scaling relation of Eq. (17) to compute it at all times for any other WD in a binary. In particular, the synchronization timescale of the lighter WD in a binary is: $${\tau }_{s,2}={\tau }_{s,1}{q}^4{\left(\frac{R_1}{R_2}\right)}^6.$$(18)

This approach provides a universal representation of tidal torques among DWDs, whose strength is captured by τ_s,ref. To explore the impact of this parameter, we consider 10 values of τ_s,ref log-uniformly distributed between 10⁻² yr and 10¹⁶ yr. The higher limit is motivated by the value of 10¹⁵ yr considered in Marsh et al. (2004), Kremer et al. (2015) and Biscoveanu et al. (2023), but, as we comment in Sect. 4, we find results to be the same for τ_s,ref ≳ 10⁸ yr. There remains a large uncertainty regarding the spin of WDs. Observations with the Kepler telescope have revealed a handful of slowly rotating WDs with rotation periods of hours to tens of hours (Hermes et al. 2017; Hernandez et al. 2024). On the other hand, a few rapidly spinning WDs with rotation periods below 40s have been observed in cataclysmic variables (Ashley et al. 2020; Lopes de Oliveira et al. 2020; Pelisoli et al. 2022), which suggests that WDs can reach high spins through accretion. To bracket uncertainties, we consider three prescriptions for the initial spins:

• ω_1,0 = ω_2,0 = 0,

• ω_1,0 = ω_2,0 = ω_o,

• ω_1,0 = 0.8ω_K,1, ω_2,0 = 0.8ω_K,2,

where ω_i,K is the break-up rate defined in Eq. (11). The first prescription represents slowly spinning WDs, the second assumes they are synchronised with the orbit since the formation of the binary, and the third prescription mimics rapidly spinning WDs (periods of ~10–100 s, in agreement with the observations of rapidly rotating WDs). Although it might not be realistic to assume that all WDs rotate rapidly at birth, exploring this scenario gives us insights into observable signatures of such type of systems. Alternatively, we could consider a scenario where only the heavier WD is highly spinning, and the lighter has virtually zero spin. It turns out that the evolution of a given binary is mostly sensitive to the spin of the heavier WD, due to its much shorter tidal synchronization timescale (see Eq. (18)). Thus, we expect the results in such a scenario would be qualitatively similar to the ones we obtain in the high-spin case, as we argue in Sect. 6.

Fig. 1 Evolution for a system with m1,0 = 0.80 M⊙ and m2,0 = 0.29 M⊙, and τs,ref = 104 yr, ω1,0 = ω2,0 = 0. The black vertical line shows the time the binary has to evolve between its formation and today. To show its long-term trajectory, we extend the evolutionary simulation beyond this timeframe.

3 Example of evolution

We computed the evolution of DWDs through numerical integration of Eqs. (2), (3), (10), and (15) with the odeint package of scipy (Jones et al. 2001). For accreting systems, the evolution can be very rapid, and the integration might fail. When this happened, we resumed the integration with a higher time resolution and kept increasing the resolution until the integration was successful. For the last stage, we decreased the resolution, in order to speed up the computation.

Typically, the separation of DWDs first shrinks due to GW radiation, and, if the binary evolves for long enough, the separation becomes sufficiently small for the lighter WD to overfill its Roche lobe, and mass transfer begins. The binary then undergoes a violent accretion episode that may cause a rapid merger (Shen 2015; Pakmor et al. 2022). If, instead, the binary survives this episode, it starts outspiralling. Following previous works (Webbink 1984; Nomoto & Iben 1985; Marsh et al. 2004; Kremer et al. 2017; Biscoveanu et al. 2023), we consider that a system becomes unstable, triggering a rapid merger, if at any point $-{\dot{m}}_d>{10}^{-2}{M}_{\odot }.{\mathrm{y}\mathrm{r}}^{-1}$.

We show in Figs. 1 and 2 the evolution of a DWD in our catalogue that undergoes stable mass transfer. For that system, τ_s,ref = 10⁴ yr and ω_1,0 = ω_2,0 = 0. The initial masses are m_a,0 = 0.80 M_⊙ and m_d,0 = 0.29 M_⊙ and the final masses m_1,f = 1.03 M_⊙ and m_2,f = 0.04 M_⊙. The vertical black line shows the time from the formation of the binary until now. In the plots shown in this section, we have evolved the systems for longer to show the long-term evolution of accreting DWDs. Fig. 1 shows the main features of the evolution of accreting systems: after the lighter WD overfills its Roche lobe, the accretion rate increases rapidly and, almost immediately, causes the distance between WDs to increase (negative ${\dot{f}}_{\text{GW}}$), reducing the orbital frequency back to a few mHz and leading to a drop-off in the accretion rate. The binary then evolves slowly in an equilibrium configuration. This equilibrium is characterised by $\dot{\mathrm{\Delta }}=0$. Deriving Eq. (1), we have $$\dot{\mathrm{\Delta }}=({\zeta }_2{R}_2-{\zeta }_{r_L}{R}_L)\frac{{\dot{m}}_2}{m_2}-R_L\frac{\dot{a}}{a},$$(19)

with ${\zeta }_2=\frac{\mathrm{d}\text{\hspace{0.17em}log}(R_2)}{\mathrm{d}\text{\hspace{0.17em}log}(m_2)}$ and ${\zeta }_{r_L}=\frac{\mathrm{d}\text{\hspace{0.17em}log}(R_L/a)}{\mathrm{d}\text{\hspace{0.17em}log}(m_2)}$. Setting $\dot{\mathrm{\Delta }}=0$ in the above equation and using Eq. (15), we get the mass transfer rate at equilibrium: $${\dot{m}}_{2,\mathrm{e}\mathrm{q}}=\frac{\frac{32}{5}\frac{G^3}{c^3}\frac{m_1{m}_{2}m}{a^4}}{1-q-\frac{m_2{j}_1}{J_{\text{orb }}}+\frac{1}{2}({\zeta }_2\left(R_2/R_L\right)-{\zeta }_{r_L})}.$$(20)

From the upper panel of Fig. 1, we can see that the equilibrium mass loss rate matches very well the actual rate. The small discrepancy is due to our neglect of sub-dominant tidal terms. We have verified that when we remove those terms from the evolution equation (Eq. (15)), the match is perfect. From the second and fourth panels of Fig. 1, we can see that most of the change in mass and frequency occurs at the very beginning of the mass transfer phase, within a few thousand years, after which they change little. The binary system loses mass only during this short time window. Subsequently, ${\dot{m}}_2$ drops by several orders of magnitude, falling below the Eddington rate, such that ${\dot{M}}_t=0$. The work of Yi et al. (2024) exploits the existence of such a simple expression for the mass-transfer rate at equilibrium, which can be written in terms of the instantaneous properties of the binary, to obtain measurements of the individual masses for accreting DWDs. In Appendix A, we discuss a case occurring in specific conditions, where the donor successively overfills and underfills its Roche lobe, leading to very pronounced oscillations in the evolution.

As shown in Fig. 2, the initial evolution is driven by GW radiation and later by mass transfer. Tidal torques initially accelerate the shrinking of the orbit. Once the WDs start exchanging mass, the heavier WD is spun up, and ω₁ > ω₀. The tidal torque on the heavier WD changes sign and helps to pull the binary apart. Thus, strong enough tidal torques can prevent binaries from getting too close and prevent them from merging rapidly.

The evolution of binaries with initially wider separations, such that they do not reach small enough separations for the donor to overflow its Roche lobe, is very similar to the early stages of the system shown in Fig. 1. The separation diminishes slowly, accompanied by a slow increase in the GW frequency, while the masses of the components remain unchanged. A noteworthy exception to this scenario occurs when WDs are born with large spins and tidal synchronization is very efficient: the binary can outspiral while far apart, without mass transfer. This happens because, for highly spinning WDs, tidal torques give a positive contribution to $\dot{a}$ (see Eq. (15)) that can be large enough to overcome GW radiation. We discuss the implications of this in the next section.

Fig. 2 Break-up of the contributions to ${\dot{f}}_{\text{GW}}$ from the different terms in Eq. (16) for the system shown in Fig. 1, zoomed in around the moment mass transfer begins. We have combined the donor and accretor contributions into a single term for the mass transfer and for the tidal torque term.

4 Results for the population of DWDs

In this section, we present the results for the whole DWD population. We excluded systems where: the accretor exceeds the Chandrasekhar mass, or $-{\dot{m}}_d>{10}^{-2}{M}_{\odot }.{\mathrm{y}\mathrm{r}}^{-1}$, or the separation goes below R₁ + R₂, as the latter is interpreted as a merger. We focus on three combinations of tidal synchronization timescales and initial spins among the thirty considered in this work, and give them the following labels:

1. ‘low spins, strong tides’: ω_1,0 = ω_2,0 = 0, τ_s,ref = 1 yr,

2. ‘low spins, weak tides’: ω_1,0 = ω_2,0 = 0, τ_s,ref = 10¹² yr,

3. ‘high spins, strong tides’: ω_1,0 = 0.8 ω_K,1, ω_2,0 = 0,8 ω_K,2, τ_s,ref = 1 yr.

Moreover, we denote by ‘GWs only’ the case where DWDs evolve solely through GW radiation. For the sake of limiting computational power, we restrict ourselves to binaries with f_GW,f above 0.1 mHz in the ‘GWs only’ case.

In Fig. 3, we show the present-day derivative of the GW frequency as a function of the final frequency itself in the ‘low spins, strong tides’ (in blue) and ‘low spins, weak tides’ (in orange) scenarios. In the upper and right panel we also plot the histograms for the ‘GWs only’ case. To make the figure more readable, we plot $F({\dot{f}}_{\mathrm{G}\mathrm{W},f})=-\mathrm{sgn}({\dot{f}}_{\mathrm{G}\mathrm{W},f}){\text{\hspace{0.17em}log}}_{10}(\mid {\dot{f}}_{\mathrm{G}\mathrm{W},f}/{10}^4{\mathrm{y}\mathrm{r}}^{-2}\mid {)}^{-1}$ rather than ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$. The minus sign and the 10⁴ yr normalisation are so that systems with positive (negative) ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ are shown as positive (negative) for the other quantity, and the inversion is so that it has the same monotonic behaviour as ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$. We show with dotted lines the contribution of accreting systems to the total histograms. Those systems typically have ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ < 0. A handful of those, at most a few tens, might have positive ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ because they just started the accretion process, and did not have time to start outspiralling, or because they are undergoing oscillations (see Appendix A).

In the f_GW,f histogram (upper panel), we observe a bump due to the accumulation of outspiralling accreting binaries around a few mHz, with a peak at ~1.3 mHz. This is not a natural limit of the evolution, but rather a typical value around which DWDs accumulate given the time they have to evolve. As described in Sect. 3, the late evolution of systems that survive the violent accretion episode at the onset of mass transfer is very slow, leading to the stacking of systems at a few mHz. Strong tidal torques help stabilise the evolution and prevent rapid mergers. Therefore, we observe that, although the bump is also observable for large values of τ_s,ref, it is more prominent for shorter tidal synchronization timescales: we find 10 times more systems at 1.3 mHz relative to the ‘GWs only’ scenario in the ‘low spins, strong tides’ case and 3 times more in the ‘low spins, weak tides’ case. The lack of events at high frequency compared to the ‘GWs only’ scenario is because binaries that reach high frequencies either merge rapidly or outspiral. Therefore, this feature is very sensitive to our criterion for the stability of mass transfer, and it is possible that this is a failure of our model. Finally, we note that the low frequency portion of the population (below 1 mHz) is pushed to higher frequencies due to tidal torques that accelerate the evolution of the systems (positive contribution to ${\dot{f}}_{\text{GW}}$, see Fig. 2), with a more pronounced effect in the ‘low spins, strong tides’ case.

Next, we compare the ‘low spins, strong tides’ (in blue) and ‘high spins, strong tides’ (in red) scenarios on Fig. 4. In the latter scenario, the initial spins of the WDs are much larger than the orbital frequency, and tidal torques give a large positive contribution to $\dot{a}$ even without mass transfer (see Eq. (15)). Therefore, tidal torques tend to slow down the evolution, and in some cases they can even dominate over GW radiation, leading the binary to outspiral. This is why, in that case, some systems are below 0.1 mHz, although the original population we evolve contains only binaries that have f_GW,f > 0.1 mHz in the ‘GWs only’ case. Some of the binaries below 0.1 mHz have positive frequency derivative because GW radiation ended up taking over again, causing them to inspiral. If a handful of DWDs were detected at low-frequency (≲0.1 mHz) with negative ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$, it would suggest that those are highly-spinning WDs and that tidal synchronization is efficient within the binary. However, as we show in Sect. 5, it might not be trivial to detect confidently such a subpopulation given the small values of |${\dot{f}}_{\mathrm{G}\mathrm{W},f}$| they have and the low signal-to-noise ratio (S/N) of sources at low frequency. We stress that for this subpopulation to form, it is not necessary that all WDs have large spin. Finally, the peak around a few mHz is less prominent in the high-spins scenario than in the low-spins one because of the slow-down of the evolution due to tidal torques: some systems do not have time to start accreting, whereas they do in the low-spins scenario.

As τ_s,ref increases, the subpopulation of low-frequency binaries with negative ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ disappears, and the initial spin has no impact. This can be seen in Fig. 5, where we show the number of accreting systems for each configuration, which is closely related to the prominence of the peak around a few mHz. For τ_s,ref ≳ 10⁸ yr, the initial spin frequency has little importance because tidal torques become negligible. For smaller tidal synchronization timescales, the number of accreting systems increases (decreases) for slowly spinning (highly spinning) WDs. When ω_1,0 = ω_2,0 = ω_o, tidal torques are initially zero, but as the separation decreases, the orbital frequency increases and the WD spins have to ‘catch up’ with it. Therefore, they tend to be smaller than ω_o, and, before the onset of mass transfer, tidal torques give a negative contribution to $\dot{a}$ (see Eq. (15)), as in the case ω_1,0 = ω_2,0 = 0. However, the absolute value of the contribution of tidal torques is less significant (because ω_o − ω_i is smaller), and some systems do not have time to start accreting within the evolution time, explaining the slight decrease in the number of accreting systems. This is a small effect, and our results support that the zero spin scenario is representative of the case where WDs are formed with small spins. For weak tidal coupling, τ_s,ref ≳ 10⁸ yr, we find ~10⁶ systems to be accreting, which represents a little more than 1% of all DWDs in the catalogue. As the tidal coupling becomes stronger, the number of accreting systems increases to ~4 × 10⁶, i.e. ~5% of all systems, for low spins, and goes down to ~3 × 10⁵, i.e. 0.4% of all systems, for high spins.

Fig. 3 Comparison of the populations in the $(f_{\mathrm{G}\mathrm{W},f},F({\dot{f}}_{\mathrm{G}\mathrm{W},f}))$ plane for the ‘low spins, strong tides’ (in blue) and ‘low spins, weak tides’ (in orange) scenarios, where $F({\dot{f}}_{\mathrm{G}\mathrm{W},f}))=-\mathrm{sgn}({\dot{f}}_{\mathrm{G}\mathrm{W},f}){\text{\hspace{0.17em}log}}_{10}(\mid {\dot{f}}_{\mathrm{G}\mathrm{W},f}/{10}^4{\mathrm{y}\mathrm{r}}^{-2}\mid {)}^{-1}$. This quantity has the same sign and monotonic behaviour as ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$. The black lines in the upper and right panels show the result in the ‘GWs only’ scenario. Accreting systems are typically those with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ < 0. We recall that we keep only DWDs that have fGW,f ≥ 0.1 mHz in the ‘GWs only’ case, which explains why the distribution goes to 0 below this limit and for small positive values of ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$.

5 Data analysis and predictions for LISA

Next, estimated the implications of the different scenarios discussed above for LISA observations. The GW signal emitted by DWDs is almost monochromatic, and takes the simple form $$h_{+}=A_0\frac{1}{2}(1+{\cos }^2(\iota ))\cos ({\varphi }_0+2\pi f_{\mathrm{G}\mathrm{W},f}t+\pi {\dot{f}}_{\mathrm{G}\mathrm{W},f}{t}^2),$$(21) $$h_{\times }=A_0\cos (\iota )\sin ({\varphi }_0+2\pi f_{\mathrm{G}\mathrm{W},f}t+\pi {\dot{f}}_{\mathrm{G}\mathrm{W},f}{t}^2),$$(22)

where $A_0=\frac{2{\mathcal{M}}_c{G}^{5/3}}{c^3{D}_L}{\left(\pi {\mathcal{M}}_c{f}_{\mathrm{G}\mathrm{W},f}\right)}^{2/3}$ gives the amplitude of the signal, ${\mathcal{M}}_c=m_1^3{m}_2^3/M_t$ is the chirp mass of the binary, D_L the luminosity distance, ι the inclination of the binary with respect to the line of sight and ϕ₀ the initial phase. We computed LISA’s response to such GW signals following Cornish & Littenberg (2007) to generate mock data, and the LISA data challenge software (LISA Data Challenge working group 2022). The sky location, which enters the LISA response, and the distance to the source were provided in the catalogue. We drew the inclination angle and the initial phase assuming isotropic orientation. We took a flat distribution for the polarisation angle, which is a quantity that enters the LISA response. We assumed a 4 yr LISA mission and took the SciRDv1 curve (LISA Science Study Team 2018) for the instrumental noise.

The detection and characterisation of individual sources in the LISA data is difficult (Littenberg & Cornish 2023; Strub et al. 2023; Zhang et al. 2021). The challenge can be attributed mostly to the very high number of different types of sources, which is expected to cause a high level of overlap between their measured GW signatures in the LISA band. For that reason, there is significant effort focusing on the developing Global Fit data analysis pipelines (Littenberg & Cornish 2023), which are based in algorithms that can simultaneously infer the parameters of each binary, as well as the actual number of binary signatures present in the data streams. As expected, these analyses are complex and often time and energy consuming. Thus, in order to characterise the number of resolvable binaries and the remaining confusion noise, we use a different approach.

Our analysis is based on the methodology presented in Karnesis et al. (2021), Timpano et al. (2006), Crowder & Cornish (2007), Nissanke et al. (2012) and Lamberts et al. (2019), which is based on an iterative procedure that classifies the sources as resolvable based on S/N criteria. The process begins (k = 0) with an initial estimate of the overall noise power spectral density S_n,k, which is the sum of the instrumental component and the confusion noise signal generated by the given population of signals. The initial estimate of the noise can be performed using different recipes, for example via a running median on the power spectral density of the data. Then, we iterate over the catalogue of sources, and classify the ones with S/N above a given threshold as resolvable. For our analysis, we chose this threshold to be equal to S/N_thres = 7. The resolvable sources are subtracted from the data, and a new iteration begins by re-estimating the overall noise S_n,k+1. The algorithm stops when no more sources can be subtracted, which can either mean that all sources in the catalogue are resolved, or that all the remaining sources are below the given threshold.

This algorithm is useful because it yields results quickly, without requiring a complete Bayesian trans-dimensional analysis. At the same time, it employs different simplifying assumptions, such as perfect data sets (no noise transients of missing data), minimal noise correlations, and perfect residuals (sources are subtracted using their true properties, which ignores uncertainties and possible correlations between the sources). Considering the above, this procedure can be seen as an optimistic upper limit to the capabilities of a full Bayesian analysis on the LISA data (Karnesis et al. 2021).

The outputs of this algorithm are summarised in Fig. 6. On the left panel, we plot the characteristic strain of the total noise power spectral density $(\sqrt{f{S}_n(f)})$, including both instrumental noise and the DWD confusion noise, in different scenarios, and on the right panel we show the number of resolvable binaries in each scenario. The hatched part of the histograms indicates the contribution of accreting binaries to all resolvable binaries. Overall, pronounced effects such as strong tides or high spins manifest observable signatures. In the ‘low spins, strong tides’ case, the confusion noise is significantly different from the ‘GWs only’ case, and in particular is much higher between 1 and 3 mHz due to accreting systems. At 1.3 mHz (the peak of accreting events), ~10⁻²% of accreting systems are resolvable, whereas ~2% of non-accreting ones at the same frequency are resolvable. This occurs because mass transfer decreases the mass ratio, and thus also the chirp mass, which reduces the GW amplitude (see Eq. (22)). As the frequency increases, this fraction increases, and above 3 mHz, accreting systems are as likely to be resolved as non-accreting ones: almost all DWDs above that frequency are resolvable. Thus, accreting systems affect mostly the confusion noise between 1 and 3 mHz. In total, ~3% of accreting binaries are resolvable, whereas ~20% of non-accreting binaries with f_GW,f > 1.3 mHz are resolvable. In the ‘low spins, strong tides’ scenario, we also predict a larger number of resolvable systems, because tidal torques accelerate the evolution of binaries, increasing the frequency at which LISA observes non-accreting systems, and therefore their GW amplitude (see Eq. (22)).

The ‘low spins, weak tides’ case yields results very similar to the ‘GWs only’ case. The bump of systems around a few mHz visible in Fig. 3 contributes to a slight increase of the confusion noise in that frequency band, and the number of resolvable systems is slightly smaller than in ‘GWs only’ scenario. The reason is that tidal effects are less effective at increasing the GW frequency and accreting systems mostly contribute to the confusion noise, but much less than in the ‘low spins, weak tides’ case, due to the smaller number of systems. The ‘high spins, strong tides’, in turn, predicts a much weaker confusion noise: the total noise differs only slightly from the instrumental noise. Moreover, it yields a much lower number of resolvable binaries: fewer sources are above 1 mHz, because of the slow-down due to tidal torques, and so they have lower S/N than in the previous scenarios. We stress that the assumption that WDs are rapidly spinning at birth is most likely not astrophysically realistic at the level of the whole population, so this result should be taken with caution. More relevant in that scenario, none of the low-frequency systems with negative ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ are resolvable, leaving for future detectors the possibility of identifying this putative subpopulation.

Several science cases with DWDs rely on measuring f_GW,f and ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ for individual systems. In order to assess how many events could be used for such studies, we use the Fisher matrix (Takahashi & Seto 2002; Vallisneri 2008) to estimate the measurement error $\mathrm{\Delta }{\dot{f}}_{\mathrm{G}\mathrm{W},f}$ (defined here as the standard deviation predicted by the Fisher matrix) on the GW-frequency derivative for resolvable systems. In Fig. 7, we show the fractional error on ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$, splitting between accreting (left panel) and nonaccreting (right panel) systems. We adopt the simple criterion $\mid \mathrm{\Delta }{\dot{f}}_{\mathrm{G}\mathrm{W},f}/{\dot{f}}_{\mathrm{G}\mathrm{W},f}\mid <1$ to decide on the measurability of ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ for a given system. In the left panel, we observe that approximately half of the resolvable systems that are not accreting have measurable GW-frequency derivative, ~5000 systems for the cases with low spins and ~1000 for the cases with high spin. However, proposed ways to measure the chirp mass and to perform tests of GR for such events rely on the evolution of DWDs being mainly driven by GWs. Even for non-accreting systems, ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ might receive significant non-GW contributions, as we now argue. Following the notation introduced in Breivik et al. (2018) and Kremer et al. (2017), we split the GW-frequency derivative into a purely gravitational (GR) and an astrophysical contribution as $${\dot{f}}_{\mathrm{G}\mathrm{W}}^{\mathrm{G}\mathrm{R}}=\frac{96}{5}\frac{f_{\mathrm{G}\mathrm{W}}^2}{\pi c^5}{\left(\pi G{\mathcal{M}}_c{f}_{\mathrm{G}\mathrm{W}}\right)}^{5/3},$$(23) $${\dot{f}}_{\mathrm{G}\mathrm{W}}^{\text{astro}}={\dot{f}}_{\mathrm{G}\mathrm{W}}-{\dot{f}}_{\mathrm{G}\mathrm{W}}^{\mathrm{G}\mathrm{R}}.$$(24)

We might interpret ${\dot{f}}_{\mathrm{G}\mathrm{W}}^{\text{astro }}$ as a systematic error, and consider it as relevant if it is larger than the statistical error. In the middle panel of Fig. 7, we show the ratio between the statistical error and ${\dot{f}}_{\mathrm{G}\mathrm{W},f}^{\text{astro }}$, restricting ourselves to systems with measurable positive GW-frequency derivative. In all three scenarios, for most of the considered systems $|\mathrm{\Delta }{\dot{f}}_{\mathrm{G}\mathrm{W},f}/f_{\mathrm{G}\mathrm{W},f}^{\text{astro }}|<1$, indicating that neglecting astrophysical effects would significantly bias the estimation of chirp mass. This occurs because systems with measurable GW-frequency derivative are in the high-frequency end of the population, where tidal torques are stronger (see Fig. 2). Within our model, in the ‘low spins, strong tides’ and ‘low spins, weak tides’ cases, we necessarily have ${\dot{f}}_{\mathrm{G}\mathrm{W}}^{\text{astro }}>0$ for the systems considered here, and the chirp mass measurement would be biased to larger values (see Eq. (23)), see also Fiacco et al. (2024). In the ‘high spins, strong tides’ case, we would instead underestimate the chirp mass. Regarding accreting systems, we find that the GW-frequency derivative is measurable for approximately half of the systems, ~5000 for the cases with strong tides and ~1000 for the cases with weak tides. Such systems could be used for the investigations suggested in Breivik et al. (2018) and Yi et al. (2024). We stress that those approaches rely, at least to some extent, on the modelling of astrophysical effects in DWDs, which therefore introduces a possible source of systematic error. On the other hand, similar approaches, for example establishing universal relations, could be used to exploit the measurement of resolvable non-accreting systems. Finally, the Fisher matrix approximation is not expected to hold exactly for all the systems in the catalogue, possibly yielding significant differences with full Bayesian analyses, but we expect our conclusions at the level of the population to hold.

Fig. 4 Same as Fig. 3, but for the ‘low spins, strong tides’ (in blue) and ‘high spins, strong tides’ (in red) scenarios. For WDs formed with high spins, tidal torques can cause the binary to outspiral even before the onset of mass transfer, explaining why some systems in red have fGW,f < 0.1 mHz.

Fig. 5 Number of accreting systems across all configurations under consideration. The left y-axis denotes the absolute count, while the right y-axis indicates the corresponding fraction relative to the initial catalogue. The DWD catalogue used as a starting point contains 7740598 systems.

Fig. 6 LISA results obtained with the subtraction algorithm in different scenarios, indicated by colours. Left: characteristic strain of the total noise power spectral density $(\sqrt{f{S}_n(f)})$, including the instrumental and the DWD confusion noise, compared to the instrumental noise alone, in grey dashed lines. Right: number of resolvable sources in each of the considered scenarios. Hatches indicate how many of those are accreting. We find 56.4%, 20.1%, and 50.7%; of the resolvable systems to be accreting in the ‘low spins, strong tides’, ‘low spins, weak tides’, and ‘high spins, strong tides’ scenarios, respectively, as indicated by the boxed numbers.

Fig. 7 Results of the Fisher matrix analysis. Left: relative error on ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ computed with the Fisher matrix for resolvable systems with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ > 0. Middle: ratio of the measurement error to the astrophysical contribution to the total GW-frequency derivative, $f_{\mathrm{G}\mathrm{W},f}^{\text{astro }}$, restricted to systems for which ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ > 0 and $\mathrm{\Delta }{\dot{f}}_{\mathrm{G}\mathrm{W},f}/{\dot{f}}_{\mathrm{G}\mathrm{W},f}<1$. Right: same as the left panel, but for resolvable systems with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ < 0. We consider quantities with fractional uncertainties less than 1 (i.e. to the left of the vertical black dashed line) to be measurable.

6 Model systematics

Given the amount of uncertainty regarding the evolution of DWDs, we test the robustness of our conclusions by considering different hypotheses in the modelling of DWDs evolution. We do not explore here the impact of modifications to the initial DWD population. Uncertainties in the star formation rate, stellar evolution and binary interactions can strongly impact properties of the population in terms of global numbers, orbital parameters, masses and positions. This will affect their later evolution and detectability by LISA (Amaro-Seoane et al. 2023).

First, we consider the approach of Gokhale et al. (2007). It is very similar to our fiducial model presented in Sect. 2, except that the specific angular momentum of the donor is not 0, but $j_2=R_2{\omega }_2^2$, and the mass lost by the system is assumed to have the same specific angular momentum as the accretor, so that ${\dot{J}}_{\text{loss }}={\dot{M}}_t{j}_1$. We find very similar results to our fiducial case. Slightly fewer systems survive accretion, due to less angular momentum being lost by the binary (smaller ${\dot{J}}_{\text{loss }}$), but the main difference is that the independency on the spins for large tidal synchronization timescales, observed on Fig. 5, no longer holds, though the difference between different initial spin scenarios is still small. This happens because the loss of angular momentum by the donor helps stabilise the evolution. Mathematically, the term $\frac{{\dot{m}}_2}{m_2}$ in Eq. (15) is now replaced by $\frac{{\dot{m}}_2}{m_2}(\frac{m_2{j}_2}{J_{\text{orb}}}+1)$, and thus it gives a more important positive contribution to $\dot{a}$, pushing the binary away, and preventing it from becoming unstable. Overall, our results confirm that the spin of the lighter WD and the model used for its specific angular momentum play a small role, unless j₂ grows with the spin and the lighter WD is rapidly spinning. Astrophysically, it is expected that the heavier WD could be rapidly spinning, due to the accretion of from the companion star before it forms the lighter WD, but it is less clear how the latter would be spun-up. The difference in the scenario where $j_2=R_2{\omega }_2^2$ and the lighter WD is rapidly spinning is however small, with ~𝒪(10³) more systems surviving accretion in the case of weak tides as compared to our fiducial model, and does not affect the main findings of the paper.

Next, given the wide uncertainty on the mass retention efficiency of WDs (Yaron et al. 2005; Bours et al. 2013; Wolf et al. 2013; Kato et al. 2017; Hillman et al. 2020), we revisit the hypothesis that WDs can accrete up to Eddington rate, and limit the maximum accretion rate of WDs to 0.1 and 0.01 of its Eddington rate. In the first case, we obtain results practically identical to the fiducial scenario. In the second case, the number of systems that survive accretion is slightly higher. This comes because the amount of mass and angular momentum lost by the systems increases, helping stabilise the evolution. The global picture remains the same: once the donor overfills its Roche lobe, a violent mass-transfer episode takes place, causing the binary either to merge or to be pushed away, the mass-transfer rate decreases dramatically, and the evolution proceeds in the equilibrium regime discussed in Sect. 3. The stacking of systems around a few mHz is a consequence of the violent episode at the onset of mass transfer and the slow subsequent evolution, and is therefore a robust prediction independent of the details of the model, provided DWDs can survive the violent mass-transfer episode.

Indeed, it has been suggested that all mass-transferring DWDs might merge (Shen 2015; Pakmor et al. 2022). To consider the impact of this hypothesis, we remove all accreting systems from our catalogue, and evaluate the stochastic signal and the resolvable binaries for this population. In this scenario, we find no strong dependency on the efficiency of tidal synchronization. This is because, once we remove all accreting systems, without distinguishing between those predicted to be stable and those predicted to be unstable, the remaining population is very similar for all values of τ_s,ref. For high-spin systems, the stochastic signal we obtain is similar to the one shown in Fig. 6 for the ‘high-spins, strong tides’ case, and the number of resolvable systems is approximately the same as the number of resolvable non-accreting systems in the ‘high-spins, strong tides’ case, given by the non-hatched part of the histogram. We stress that this latter result is not as trivial as it might seem at first glance. When changing the population, the noise due to the stochastic background changes, changing the criterion within our algorithm to decide if a source is resolvable or not. For the two other scenarios, the number of resolvable systems is similar to the nonaccreting portion of the ‘low-spins, weak tides’ and the stochastic signal is very similar to the one in that scenario, with a slight decrease around a few mHz, due to the absence of accreting systems.

7 Conclusions

In order to prepare for the analysis and astrophysical interpretation of LISA data, it is crucial to bracket the astrophysical predictions for the population of DWDs LISA will observe. In this study, we investigated the influence of tidal torques and mass transfer on the DWD population and their implications for the confusion noise DWDs generate and the number of resolvable binaries.

Starting from a mock catalogue of Galactic DWDs at their formation, predicted by a population synthesis code, we evolved these binaries to the present day using a semi-analytical model that includes the leading order effect of GW radiation, but also tidal torques and mass transfer within the binaries, as well as angular momentum loss due to non-conservation of mass when the accretion rate becomes too large. For binaries undergoing mass transfer, we applied the same simplified criterion found in the literature to determine their stability, which is based on a limit value for the rate of mass lost by the donor. Systems failing this criterion were assumed to merge rapidly and are removed from the catalogue, while those passing it are observed to outspiral and accumulate around a few mHz.

In this model, tidal torques depend on the spins of the WDs and on the tidal synchronization timescale. We define the value of the latter for a reference system at the moment it overfills its Roche lobe, and employ scaling relations based on the mass ratio and separation of the binary to compute it for all systems at any time, facilitating a universal representation of tidal torques. We performed simulations using ten different values of this parameter and explored the uncertainty regarding the spin of WDs at birth by considering three scenarios: one where WDs are born with zero-spin, one where they are already synchronised with the orbit, and one where they are rapidly spinning. We find that the first two cases yield similar results. After generating the mock population, we split the contribution of GWs emitted by DWDs in LISA into a stochastic signal, characterised by its power spectral density, and resolvable binaries.

Our findings indicate that weak tidal torques (long timescales) have minimal impact on the population, and incorporating tidal torques and mass transfer into the evolutionary model results in only marginal adjustments compared to scenarios driven solely by GW radiation. Conversely, strong tidal torques can stabilise the evolution and prevent rapid mergers. For WDs born with low spins, we observe a significant accumulation of binaries around a few mHz, which substantially impacts the confusion noise in that frequency region compared to scenarios with GWs alone. Additionally, there is a higher number of resolvable systems in this case. Conversely, for WDs born with large spins, tidal torques tend to slow down the evolution and can even overcome GW radiation, causing the binaries to outspiral even before mass transfer begins. Consequently, the confusion noise due to DWDs becomes almost negligible and very few binaries are resolvable. It is important to note that this latter scenario is likely not representative of the entire population, and predictions for LISA based on it should be interpreted cautiously. On the other hand, observing binaries at very low frequencies (≲ 0.1 mHz) with negative frequency derivatives could indicate the presence of a subpopulation of DWDs born with large spins and strong tidal torques. However, due to the low S/N of low-frequency DWDs, we estimate that LISA may not have the sensitivity to detect such a subpopulation, leaving this possibility for future detectors. In a nutshell, strong effects (i.e. strong tides and/or high spins) yield observable signatures.

Using the Fisher matrix approximation, we estimated the error on the GW-frequency derivative for resolvable systems. Although the latter is limited, in principle, to the high S/N regime (Vallisneri 2008), it allows us to estimate the parameters errors of each binary without having to perform a full trans-dimensional Bayesian analysis. We find ~5000 and ~1000 accreting systems to have measurable negative GW-frequency derivative in the case of strong and weak tides, respectively By combining with assumptions on the evolution of DWDs, it should be possible to measure the mass of such systems (Breivik et al. 2018; Yi et al. 2024). Regarding non-accreting systems, we find ~5000 and ~1000 of them to have measurable positive GW-frequency derivative in the case of low and high natal spins, respectively. However, for all of them, the systematic error resulting from neglecting astrophysical effects, due to tidal torques in our model, is larger than the statistical error. This renders them unsuitable for chirp mass measurements and tests of GR based on the joint measurement of the gravitational wave frequency and its first derivative in DWDs, unless astrophysical information is incorporated into these measurements, for example by means of universal relations.

We find our results to be robust to different assumptions regarding the evolution of DWDs. The main uncertainty resides in the stability of mass-transferring binaries. in our study, systems were labelled as unstable under mass transfer when the rate of mass lost by the donor exceeded 10⁻²M_⊙.yr⁻¹. We expect our results to change monotonically with this mass-transfer rate limit: a larger limit rate would lead to more stable accreting binaries and vice versa. We encountered numerical challenges in solving for the evolution within this high-rate regime even within our semi-analytical model, suggesting this is indeed an uncertain regime. However, given the importance of mass transfer on the population of DWDs, a deeper understanding of the stability conditions is essential, in the lines of Chen et al. (2022) The lack of systems at high frequencies (f_GW,f > 20mHz) is also a consequence of our simplistic treatment: we removed all systems that exceed the limit mass-transfer rate from the catalogue. A better categorisation and treatment of unstable systems could modify this feature. Additionally, for the high-frequency DWDs, which are the ones with measurable frequency derivative, non-linearities in tidal torques should be accounted for (Yu et al. 2020). Moreover, we restricted ourselves to a single population of DWDs at formation predicted by our population synthesis code. However, there are also numerous uncertainties regarding star formation rate, stellar evolution and binary interactions that would affect the observed population of DWDs. Finally, we did not account for differences in the composition of WDs. It would be interesting to assess how modifications to the mass-radius relation (Hamada & Salpeter 1961; Suh & Mathews 2000; Pane et al. 2000; Karinkuzhi et al. 2024) would impact the observed population, creating perhaps distinctive signatures that would allow identifying the type of WDs within a binary. We leave these endeavours for future research.

Our results complement the findings of Scaringi et al. (2023), where the authors found that cataclysmic variables, formed of a WD accreting from a main-sequence or sub-stellar star, lead to measurable deviations in the confusion noise below 1 mHz Together, these studies highlight the need for flexible tools in LISA data analysis, in order to account for all possibilities in the widely unknown population of compact objects in the Galaxy. Our exploration provides insights into potential parametrisations of the DWD population for population analysis with LISA. We suggest representing the population of resolvable DWDs as a power-law contribution in frequency with positive ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ and a Gaussian around a few mHz with negative ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$. For the stochastic signal, a flexible function accommodating the various possibilities demonstrated in our study would be suitable, possibly employing the linear interpolation scheme proposed with variable number of knots proposed in Toubiana et al. (2023). In order to gain informative power, it would be valuable to establish relationships between the properties of resolvable and non-resolvable populations, a task we defer to future investigations.

Acknowledgements

We are thankful to V. Korol for her valuable suggestions during the preparation of this work. We are also thankful to the efforts of the LDC Working Group for the software availability and help. N.Karnesis acknowledges the funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101065596. Astrid Lamberts is supported by the ANR COSMERGE project, grant ANR-20-CE31-001 of the French Agence Nationale de la Recherche. This work was supported by the ‘Programme National des Hautes Energies’ (PNHE) of CNRS/INSU co-funded by CEA and CNES.

Appendix A Oscillations in the evolution

For certain combinations of WD masses and tidal synchronization timescale, we observe very pronounced oscillations in the binary separation due to the donor successively overfilling and ‘underfilling’ its Roche lobe. A similar behaviour is discussed in Marsh et al. (2004); Gokhale et al. (2007); Kremer et al. (2015). This oscillatory behaviour is illustrated in Fig. A.1, where we show the evolution of ${\dot{f}}_{\mathrm{G}\mathrm{W}}$ for the same system as in Fig. 2, but taking a shorter tidal synchronization timescale, τ_s,ref = 10² yr. The explanation behind these oscillations is as follows. After the donor first overfills its Roche lobe and mass transfer starts, the distance between WDs increases, at a larger rate than in the previous case, due to the larger negative (positive for the separation) contribution of tidal terms (see upper panel of Fig. A.1), the Roche lobe expands faster than the donor, which no longer overfills its Roche lobe. GWs take over and bring the system closer again, until the donor overfills its Roche lobe and mass transfer resumes. This succession of events keeps on repeating until the donor is light enough and its radius expands faster than the Roche lobe, so that mass transfer is not interrupted. In total, this oscillatory stage would last a few tens of Gyr, alternating oustpiralling regimes of ~ 1 Gyr and inspiralling regimes of ~ 5 Gyr.

Fig. A.1 Same as Fig. 2, but with τs,ref = 102 yr, and we show in parallel the evolution of the Roche lobe overflow.

If we consider an even shorter τ_s,ref, we no longer observe these oscillations. The reason for this is as follows. Before the onset of mass transfer, tidal torques are very efficient at synchronising the WDs with the orbit, as can be seen in Fig. A.2. Then, mass transfer starts, leading to an additional torque (first term in Eq. (10)), which desynchronises the WDs (in this example, mostly the accretor). The tidal torque then becomes very strong, due to the small synchronization timescale, and drives the binary apart. Once the tidal term becomes sub-dominant, the donor is light enough so that the Roche lobe does not expand faster than its radius, preventing the oscillatory behaviour described above.

Fig. A.2 Same as Fig 2, but with τs,ref = 1 yr, and we show in parallel the evolution of the accretor and the donor spins, as well as the orbital frequency.

References

All Figures

	Fig. 1 Evolution for a system with m1,0 = 0.80 M⊙ and m2,0 = 0.29 M⊙, and τs,ref = 104 yr, ω1,0 = ω2,0 = 0. The black vertical line shows the time the binary has to evolve between its formation and today. To show its long-term trajectory, we extend the evolutionary simulation beyond this timeframe.
In the text

	Fig. 2 Break-up of the contributions to ${\dot{f}}_{\text{GW}}$ from the different terms in Eq. (16) for the system shown in Fig. 1, zoomed in around the moment mass transfer begins. We have combined the donor and accretor contributions into a single term for the mass transfer and for the tidal torque term.
In the text

Fig. 3 Comparison of the populations in the $(f_{\mathrm{G}\mathrm{W},f},F({\dot{f}}_{\mathrm{G}\mathrm{W},f}))$ plane for the ‘low spins, strong tides’ (in blue) and ‘low spins, weak tides’ (in orange) scenarios, where $F({\dot{f}}_{\mathrm{G}\mathrm{W},f}))=-\mathrm{sgn}({\dot{f}}_{\mathrm{G}\mathrm{W},f}){\text{\hspace{0.17em}log}}_{10}(\mid {\dot{f}}_{\mathrm{G}\mathrm{W},f}/{10}^4{\mathrm{y}\mathrm{r}}^{-2}\mid {)}^{-1}$. This quantity has the same sign and monotonic behaviour as ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$. The black lines in the upper and right panels show the result in the ‘GWs only’ scenario. Accreting systems are typically those with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ < 0. We recall that we keep only DWDs that have fGW,f ≥ 0.1 mHz in the ‘GWs only’ case, which explains why the distribution goes to 0 below this limit and for small positive values of ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$.

In the text

	Fig. 4 Same as Fig. 3, but for the ‘low spins, strong tides’ (in blue) and ‘high spins, strong tides’ (in red) scenarios. For WDs formed with high spins, tidal torques can cause the binary to outspiral even before the onset of mass transfer, explaining why some systems in red have fGW,f < 0.1 mHz.
In the text

	Fig. 5 Number of accreting systems across all configurations under consideration. The left y-axis denotes the absolute count, while the right y-axis indicates the corresponding fraction relative to the initial catalogue. The DWD catalogue used as a starting point contains 7740598 systems.
In the text

Fig. 6 LISA results obtained with the subtraction algorithm in different scenarios, indicated by colours. Left: characteristic strain of the total noise power spectral density $(\sqrt{f{S}_n(f)})$, including the instrumental and the DWD confusion noise, compared to the instrumental noise alone, in grey dashed lines. Right: number of resolvable sources in each of the considered scenarios. Hatches indicate how many of those are accreting. We find 56.4%, 20.1%, and 50.7%; of the resolvable systems to be accreting in the ‘low spins, strong tides’, ‘low spins, weak tides’, and ‘high spins, strong tides’ scenarios, respectively, as indicated by the boxed numbers.

In the text

Fig. 7 Results of the Fisher matrix analysis. Left: relative error on ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ computed with the Fisher matrix for resolvable systems with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ > 0. Middle: ratio of the measurement error to the astrophysical contribution to the total GW-frequency derivative, $f_{\mathrm{G}\mathrm{W},f}^{\text{astro }}$, restricted to systems for which ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ > 0 and $\mathrm{\Delta }{\dot{f}}_{\mathrm{G}\mathrm{W},f}/{\dot{f}}_{\mathrm{G}\mathrm{W},f}<1$. Right: same as the left panel, but for resolvable systems with ${\dot{f}}_{\mathrm{G}\mathrm{W},f}$ < 0. We consider quantities with fractional uncertainties less than 1 (i.e. to the left of the vertical black dashed line) to be measurable.

In the text

	Fig. A.1 Same as Fig. 2, but with τs,ref = 102 yr, and we show in parallel the evolution of the Roche lobe overflow.
In the text

	Fig. A.2 Same as Fig 2, but with τs,ref = 1 yr, and we show in parallel the evolution of the accretor and the donor spins, as well as the orbital frequency.
In the text