# Applying the metallicity-dependent binary fraction to double white dwarf
formation:
Implications for LISA

Sarah Thiele Department of Physics and Astronomy, University of British
Columbia, 6224 Agricultural Road, Vancouver, BC, V6T 1Z1, Canada  Canadian
Institute for Theoretical Astrophysics, University of Toronto, 60 St. George
Street, Toronto, Ontario, M5S 1A7, Canada  Katelyn Breivik Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY,
10010, USA  Canadian Institute for Theoretical Astrophysics, University of
Toronto, 60 St. George Street, Toronto, Ontario, M5S 1A7, Canada  Robyn E.
Sanderson Department of Physics and Astronomy, University of Pennsylvania, 209
South 33rd Street, Philadelphia, PA 19104, USA  Center for Computational
Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY, 10010, USA

###### Abstract

Short-period double white dwarf (DWD) binaries will be the most prolific
source of gravitational waves (GWs) for the Laser Interferometer Space Antenna
(LISA). Not only will tens of thousands of DWDs be individually resolved, but
DWDs with GW frequencies below $\sim 1\,\rm{mHz}$ will be the dominant
contributor to a stochastic foreground caused by confusion from overlapping GW
signals, limiting the detectability of individual sources of all kinds.
Population modeling of Galactic DWDs typically assumes a standard binary
fraction of $50\%$. However, recent observations have shown that the binary
fraction of close ($P_{\text{orb}}$ $\leq 10^{4}$ days) solar-type stars
exhibits a strong anti-correlation with metallicity. In this study we perform
the first simulation of the Galactic DWD population observable by LISA which
incorporates an empirically-derived metallicity-dependent binary fraction. We
simulate DWDs using the binary population synthesis suite COSMIC and
incorporate a metallicity-dependent star formation history to create a
Galactic population of short-period DWDs. We compare two models: one which
assumes a metallicity-dependent binary fraction, and one with a binary
fraction of 50%. We find that while metallicity impacts the evolution and
intrinsic properties of our simulated DWD progenitor binaries, the LISA-
resolvable populations of the two models remain roughly indistinguishable.
However, the size of the total Galactic DWD population orbiting in the LISA
frequency band is reduced by more than half when accounting for a metallicity-
dependent binary fraction. This effect serves to lower the confusion
foreground, effectively increasing the sensitivity for detecting all types of
low-frequency LISA sources. We repeat our analysis for three different
assumptions for Roche-lobe overflow interactions and find the population
reduction to be robust when a metallicity-dependent binary fraction is
assumed.

Binary stars — Stellar evolution — GW astronomy

††software: astropy (Astropy Collaboration et al., 2013; Price-Whelan et al.,
2018); COSMIC (Breivik et al., 2020b); LEGWORK (Wagg et al., 2021); matplotlib
(Hunter, 2007); numpy (van der Walt et al., 2011); pandas (Wes McKinney, 2010;
pandas development team, 2020); scipy (Jones et al., 2001)

##  1 Introduction

Most stars in the Galaxy will end their lives as white dwarfs. Of the stars
which are born with a binary companion, many will undergo interactions which
bring the two stars closer together, eventually forming a close double white
dwarf (DWD). Close DWDs, with orbital periods shorter than $\sim
2.75\,\rm{hr}$ are the largest source by number of mHz gravitational waves
(GWs) in the Galaxy (e.g., Amaro-Seoane et al., 2017). The Laser
Interferometer Space Antenna (LISA) is expected to resolve at least $10^{4}$
individual DWD binaries in the Milky Way and will also observe GW emission
from the entire Galactic DWD population through the unresolved foreground
created by overlapping signals at sub-mHz frequencies (e.g., Nelemans et al.,
2001; Ruiter et al., 2010; Nissanke et al., 2012; Yu & Jeffery, 2013; Korol et
al., 2017; Lamberts et al., 2019; Breivik et al., 2020b). The resolved
population will enable the study of several important aspects of binary
evolution like the strength of tides (Valsecchi et al., 2012) and the
stability of mass transfer in DWD systems (e.g., Marsh et al., 2004; Shen,
2015; Gokhale et al., 2007; Sepinsky & Kalogera, 2014; Kremer et al., 2015) as
well as provide a probe of Galactic structure (Korol et al., 2019) and the
Local Group (Korol et al., 2018). The shape and strength of the Galactic DWD
foreground can also be used as a tool to study the structure of the Milky Way
(Benacquista & Holley-Bockelmann, 2006; Breivik et al., 2020a) as well as the
separation distribution of close DWDs (Korol et al., 2021).

When viewed strictly as a source of noise, the unresolved Galactic DWD
foreground is the dominant noise source for LISA in the sub-mHz part of LISA’s
frequency range. This extra noise above the detector noise floor affects the
detection of all other LISA sources including extreme mass ratio inspirals
(e.g., Berti et al., 2006; Barack & Cutler, 2007; Babak et al., 2017; Moore et
al., 2017), merging black holes with masses between $10^{4}$–$10^{7}$ M⊙
(e.g., Klein et al., 2016; Bellovary et al., 2019), and cosmological GW
backgrounds (e.g., Bartolo et al., 2016; Caprini et al., 2016; Caldwell et
al., 2019). For sources which have signals buried by the Galactic DWD
population, the foreground must be carefully analyzed and subtracted (Adams &
Cornish, 2014; Cornish, 2020; Littenberg et al., 2020; Boileau et al., 2021).
The number of resolved DWDs and the height of the unresolved DWD foreground
are a direct consequence of the number of DWD progenitors which form and
evolve over the Milky Way’s history.

While the binary fraction remains approximately constant across a large
metallicity range ($-1.5\leq$ [Fe/H] $<0.5$) for wide binaries, close OB
stars, and the stellar Initial Mass Function (IMF) (Moe & Di Stefano, 2017;
Moe et al., 2019), the binary fraction for solar-type star systems with
orbital period $P_{\text{orb}}\leq 10^{4}$ days (separation $a\leq 10$ AU)
shows a strong anti-correlation with metallicity (e.g. Badenes et al., 2018;
Moe et al., 2019; Mazzola et al., 2020; Price-Whelan et al., 2020). Because
close DWDs are the remnants of close, solar-type binary stars, this anti-
correlation plays an important role in the formation, evolution, and
characteristics of the DWD population that LISA will observe.

To date, population synthesis studies of the Galactic population of close DWDs
have either assumed a $100\%$ binary fraction or a $50\%$ binary fraction,
such that for every three stars formed, two reside in a binary system
(Nelemans et al., 2001; Yu & Jeffery, 2013; Korol et al., 2017; Lamberts et
al., 2019). In this study, we investigate the effects of a metallicity-
dependent binary fraction on the formation and evolution of DWDs. To this end,
we create synthetic present-day Milky Way-like galaxies of DWDs and
specifically select systems with GW signals that may be observable by the
space-based detector LISA. Throughout, we make comparisons between the
standard assumption of a constant $50\%$ initial binary fraction (hereafter
model F50) and one with a metallicity-dependent binary fraction (hereafter
model FZ).

In Section 2 we discuss our binary evolution assumptions used in simulating
DWD populations and detail the process to produce present-day synthetic Milky-
Way-like galaxies. In Section 3 we review the derivation of LISA detectability
for circular DWD populations at mHz frequencies. In Section 4 we detail
results showing how a metallicity-dependent binary fraction affects the
formation and evolution of DWD populations assuming a fiducial set of binary
evolution assumptions. In Section 5 we detail the metallicity dependence of
the LISA DWD foreground and resolved population. Finally, we repeat our
analysis for three sets of binary evolution assumptions to show a robust
reduction of the size of the Galactic DWD foreground in LISA for populations
simulated with a metallicity-dependent binary fraction in Section 6 and
conclude in Section 7.

##  2 Simulating a Galactic DWD population

In this section we describe the setup of our DWD simulations using the binary
population synthesis suite COSMIC, and the process to scale these simulations
to create Milky Way-like galaxies using the star formation history of galaxy
m12i in the Latte suite of the FIRE-2 simulations (Wetzel et al., 2016;
Hopkins et al., 2018) and stellar positions assigned according to the the
Ananke framework (Sanderson et al., 2020).

###  2.1 Binary population models

We simulate the evolution of DWD progenitor populations using
COSMIC111https://cosmic-popsynth.github.io, an open-source Python-based rapid
binary population synthesis suite which employs single and binary star
evolution using SSE/BSE (Hurley et al., 2000, 2002). Several modifications
have been added to COSMIC which incorporate updates for massive star evolution
and binary interactions. For a detailed description of these modifications see
Breivik et al. (2020b). COSMIC has been used in several studies to examine the
effects of binary evolution on binary populations from blue stragglers (Leiner
& Geller, 2021) and heartbeat stars (Jayasinghe et al., 2021), to white dwarf
populations (Kremer et al., 2017; Breivik et al., 2018; Kilic et al., 2021),
to merging compact object populations in isolated binaries (Zevin et al.,
2020a, b, 2021; Wong et al., 2021; Mandhai et al., 2021) and in dynamical
environments around super-massive black holes (Stephan et al., 2019; Wang et
al., 2021).

COSMIC is especially useful for efficient generation of large populations of
compact binaries. Instead of choosing a fixed number of binary stars for each
simulation, COSMIC iteratively simulates populations until parameter
distributions of the binary population converge to a stable shape as more
binaries are added. This process is quantified through the $match$ parameter
inspired by matched filtering techniques (e.g. Eq. 6 of Chatziioannou et al.,
2017) defined as

$match=\frac{\sum_{k=1}^{N}P_{k,i}P_{k,i+1}}{\sqrt{\sum_{k=1}^{N}P_{k,i}P_{k,i}\sum_{k=1}^{N}P_{k,i+1}P_{k,i+1}}},$
(1)

where $P_{k,i}$ represents the height of bin $k$ on the $i\rm{th}$ iteration
(Breivik et al., 2020b). In this study, we simulate binaries until
$\log_{10}(1-match)\leq-5$ for the masses and orbital periods of each DWD
population at the formation of the second WD. Since all DWD progenitor
binaries simulated with COSMIC are circularized through mass transfer or tides
before the second WD forms (e.g. Marsh et al., 2004; Gokhale et al., 2007;
Sepinsky & Kalogera, 2014; Kremer et al., 2015), we do not consider
convergence of DWD eccentricities.

The masses and orbital periods at the formation of the second-formed WD span a
wide range depending on the WD binary component types, thus we consider four
DWD combinations: two helium WDs (He + He), a carbon-oxygen WD orbiting a
helium WD (CO + He), two carbon-oxygen WDs (CO + CO), and an oxygen-neon
orbiting a helium, carbon-oxygen, or oxygen-neon WD (ONe + X). For each DWD
type we simulate a grid of $15$ metallicities spaced uniformly in
$\log_{10}(Z)$ between $Z=10^{-4}$ to $0.03$, to account for the limits of the
Hurley et al. (2000) stellar evolution tracks employed in COSMIC. This results
in a total of $60$ populations across all DWD types and metallicities. The
output of COSMIC contains information limited to intrinsic binary properties
like mass and orbital period. External parameters like Galactic position and
orientation are assigned in a post-processing scheme which uses metallicity-
dependent positions and ages from the Ananke framework of galaxy m12i from the
Latte Suite of the FIRE-2 simulations (see Section 2.3 for details).

We assume that the Zero Age Main Sequence (ZAMS) masses, orbital periods and
eccentricities for each binary are independently distributed. We choose
primary masses following Kroupa (2001), a flat mass ratio distribution (Mazeh
et al., 1992; Goldberg & Mazeh, 1994), a log-uniform period distribution
following Opik’s Law, and a uniform eccentricity distribution following Geller
et al. (2019). We initialize all binaries with the same evolution time of
$13.7\,\rm{Gyr}$ to capture all potential evolution within a Hubble time. We
assume a $100\%$ binary fraction in our COSMIC simulations to reduce
computation time and scale the simulations to models which assume a constant
$50\%$ binary fraction (model F50) or a metallicity-dependent binary fraction
(model FZ) in a post-processing scheme.

Since we are primarily interested in the effects of a metallicity-dependent
binary fraction, we use a single set of assumptions for binary interactions.
Our choices follow the COSMIC defaults described in Breivik et al. (2020b)
except for the treatment of Roche-lobe overflow (RLO). The stability of RLO
mass transfer is determined using critical mass ratios resulting from radius-
mass exponents (Webbink, 1985; Hurley et al., 2002), where the critical mass
ratio is defined as the ratio of the donor to accretor mass. We assume
critical mass ratios following Claeys et al. (2014) which reduce the standard
critical mass ratio assumptions from Hurley et al. (2002) for main sequence
(MS) donors by $\sim 50\%$ from $3$ to $1.6$ based on the models of de Mink et
al. (2007) and treat WD accretors separately following the models of Soberman
et al. (1997). We increase the mass loss rate from the donor following
Equation 11 of Claeys et al. (2014). The amount of mass lost during RLO from
the donor is limited by the overflow factor of the donor radius to its Roche
radius following Hurley et al. (2002). The amount of mass accepted by the
accretor is limited to $10$ times the accretor’s mass divided by the
accretor’s thermal timescale. Finally, for RLO mass loss which becomes
unstable and leads to common envelope evolution (CEE) we assume that the
donor’s binding energy is calculated according to the fits detailed in
Appendix B of Claeys et al. (2014) and that orbital energy is deposited with
$100\%$ efficiency into unbinding the common envelope ($\alpha=1$).

###  2.2 Metallicity-dependent binary fraction

We fit the results presented in Moe et al. (2019) using linear regression to
obtain a piecewise relation between the metallicity, [Fe/H] and binary
fraction $f_{\rm{b}}$ as

$f_{\rm{b}}=\begin{cases}-0.0648\cdot[\rm{Fe/H}]+0.3356,&[\rm{Fe/H}]\leq-1.0\\\
-0.1977\cdot[\rm{Fe/H}]+0.2025,&[\rm{Fe/H}]>-1.0\end{cases}$

We convert between [Fe/H] and metallicity Z, assuming all stars have solar
abundance such that

$[\rm{Fe/H}]=\log_{10}\Big{(}\frac{Z}{Z_{\odot}}\Big{)},$ (2)

where we assume $Z_{\odot}=0.02$.

###  2.3 A metallicity-dependent SFH: Convolving with the FIRE-2 models

To create Milky Way-like galaxies which integrate the metallicity-dependent
binary fraction, we use the metallicity-dependent ages and positions of galaxy
m12i from the “Latte” suite of the FIRE-2 simulations (Hopkins, 2015; Wetzel
et al., 2016; Hopkins et al., 2018) to create synthetic, Milky Way-like DWD
populations.

The m12i galaxy provides particle mass resolution of 7070 M⊙ per star
particle. Each star particle has an associated metallicity, position, and age,
which is combined with the output of COSMIC to assign DWDs to each star
particle by matching its metallicity to our COSMIC metallicity grid. The
positions of each DWD are assigned using the Ananke framework since multiple
DWD binaries can form within a single star particle. Specifically, we use an
epanechnikov kernel where the kernel size is inversely proportional to the
local density to assign the radial component of spherically symmetric offsets
from the center of each star particle following Sanderson et al. (2020).

Figure 1: The metallicity-dependent binary fraction for close solar-type
binaries with P $<10^{4}$ days (black) plotted against the logarithm of Z/Z⊙,
with scatter points denoting the location of the metallicity grid used in our
COSMIC simulations. The secondary axis shows the amount of mass in stars
formed within each metallicity bin from galaxy m12i of the “Latte” suite in
the FIRE-2 simulations as a red histogram. The amount of stellar mass formed
with super-solar metallicities dominates the distribution. Note that the
primary axis shows $f_{\rm{b}}(Z)$ in linear scale, and the secondary axis
shows amount of mass formed in log-scale. The opposing trends of these two
distributions compete throughout this study.

The metallicity-dependent binary fraction is shown in black in Figure 1 along
with the mass in star particles from galaxy m12i, shown in red, as a function
of our metallicity grid. The binary fraction, $f_{\rm{b}}$, drops drastically
across metallicity while the mass formed in m12i increases significantly.
These two opposing trends compete throughout this study along with the impact
of metallicity on single star evolution to form the final numerical
distribution of systems in our DWD populations.

Since our COSMIC simulations assume a binary fraction of $f_{\rm{b}}=1$, we
scale the amount of mass sampled at ZAMS required to produce our COSMIC-
generated population of DWDs ($M_{\text{b,ZAMS}}$) to the proper amount of
mass sampled in single and binary stars ($M_{\text{ZAMS,sim}}$) for each
binary fraction and model. We do this by sampling single stars and primary
masses of binary stars from the Kroupa (2001) IMF and sampling secondary
masses of the binary stars from a uniform mass distribution, where the number
of binaries is determined by the binary fraction model (F50 or FZ). From this
sample, we obtain the ratio of mass in single stars to the mass in binary
stars as a function of the binary fraction, $R(f_{\rm{b}})$. For model F50,
the ratio is a constant $R(f_{\rm{b}})=0.64$. For model FZ, the ratio
increases from 0.68 at low metallicities to 3.17 at solar metallicity
indicating less mass in binaries relative to single stars. The total amount of
ZAMS mass in single and binary stars is then
$M_{\text{ZAMS,sim}}=M_{\text{b,ZAMS}}(1+R(f_{\rm{b}}))$.

Once we determine the total ZAMS mass required to produce our simulated
population for a given metallicity, the number of DWDs formed per unit solar
mass at metallicity $Z_{i}$ is

$n_{\rm{DWD}}(Z_{i})=\frac{N_{\rm{DWD,sim}}(Z_{i})}{M_{\rm{ZAMS,sim}}(Z_{i})}.$
(3)

The number of DWD’s per m12i star particle at metallicity $Z_{i}$ is then

$N_{\rm{DWD},\star}(Z_{i})=n_{\rm{DWD}}(Z_{i})\,M_{\star},$ (4)

where $M_{\star}=7070\,M_{\odot}$ is the mass per m12i star particle. Since
$N_{\rm{DWD},\star}(Z_{i})$ is not an integer, we treat the decimal component
as the probability that the star particle contains an extra DWD in addition to
the integer number. For each star particle, we sample with replacement
$N_{\rm{DWD},\star}(Z_{i})$ DWDs from the corresponding simulated COSMIC
population at that metallicity and assign the ZAMS birth time of each DWD to
the formation time of the star particle. For most DWD types there is more than
one DWD binary system assigned to each m12i star particle.

If the DWD formation time is less than the age of the star particle, we evolve
the DWD over the remaining time between its formation and star particle age,
$t_{\rm{evol}}$, to produce the present-day population. Once a DWD is formed,
we assume that the binary evolves only due to the emission of GWs. Due to
tidal effects and mass transfer between their progenitor binaries, all DWDs in
our simulations are circular, thus eccentricity does not need to be
considered. The orbital evolution over the time $t_{\rm{evol}}$ is then simply
defined according to Peters (1964) as

$a_{f}=(a_{i}-4\beta t_{\text{evol}})^{1/4},$ (5)

where $a_{i}$ is the DWD separation at formation, and

$\beta=\frac{64G^{3}}{5c^{5}}M_{1}M_{2}(M_{1}+M_{2})$ (6)

is constant throughout DWD evolution (Peters, 1964).

We discard any DWDs for which the sum of their ZAMS birth time, given by the
star particle formation time, and DWD formation time is larger than the age of
the star particle since the system will not have evolved long enough to become
a DWD at present. We further discard any DWDs for which the lower mass WD
overflows it’s Roche lobe before present day, because the outcomes of these
interactions are highly uncertain and their treatment is outside the scope of
this work (e.g., Shen, 2015; Kremer et al., 2017). The separation at which the
lower mass WD overflows its Roche Lobe is defined as

$a_{\text{RLO,}\ell}=R_{\ell}\frac{0.6q_{\ell}^{2/3}+\ln{(1+q_{\ell}^{1/3})}}{0.49q_{\ell}^{2/3}}$
(7)

where $R_{\ell}$ is the radius of the lower mass WD and
$q_{\ell}=M_{\ell}/M_{h}$ is the ratio of the lower- to higher-mass WD
components. (Eggleton, 1983). We define the radius of a WD following Tout et
al. (1997); Hurley et al. (2000) as

$R_{\rm{WD}}=\max\Bigg{(}R_{\text{NS}},0.0115\sqrt{\left(\frac{M_{\text{Ch}}}{M}\right)^{2/3}-\left(\frac{M}{M_{\text{Ch}}}\right)^{2/3}}\Bigg{)}$
(8)

where $R_{\text{NS}}=1.4\cdot 10^{-5}$ R⊙ is the radius of a neutron star,
$M_{\text{Ch}}=1.44$ M⊙ is the Chandrasekhar limit for the mass of a stable
WD, and $M$ is the mass of the WD in solar masses.

For the non-discarded systems, we log the present-day separations from which
the present-day orbital frequency $f_{\text{orb}}$ can be found using Kepler’s
third law. The GW frequency is then $f_{\text{GW}}=2f_{\text{orb}}$.

##  3 LISA detectability

We use LEGWORK222https://legwork.readthedocs.io (Wagg et al., 2021) to
determine the detectability of our simulated DWD populations for sources with
GW frequencies $f_{\rm{GW}}>10^{-4}\,\rm{Hz}$. LEGWORK calculates the
position-, orientation-, and angle-averaged signal to noise ratio (SNR) for
inspiraling GW sources closely following the derivations of Flanagan & Hughes
(1998) and using the LISA noise power spectral density (PSD) of Robson et al.
(2019).

To lowest order in the post-Newtonian expansion, the frequency evolution of
circular orbits for quadrupole GW emission is defined as

$\dot{f}_{n}=\frac{48n}{5\pi}\frac{(G\mathcal{M}_{c})^{5/3}}{c^{5}}(2\pi
f_{\text{orb}})^{11/3}.$ (9)

We classify DWDs as evolving, or “chirping”, when $\dot{f}_{n}\geq
1/t_{\text{obs}}^{2}$. For evolving sources, the SNR is

$\langle\rho\rangle^{2}_{\text{circ,evol}}=\int_{f_{0}}^{f_{1}}df\frac{h_{c}^{2}}{f^{2}S_{n}(f)}$
(10)

where $h_{c}$ is the characteristic strain of the system, $S_{n}(f)$ is the
LISA sensitivity curve of Robson et al. (2019), and the frequency limits are
determined by the orbital evolution over the observation time,
$T_{\rm{obs}}=4\,\rm{yr}$. The characteristic strain for circular orbits is

$h_{c}^{2}=\frac{2^{2/3}}{3\pi^{4/3}}\frac{(G\mathcal{M}_{c})^{5/3}}{c^{3}D_{L}^{2}}\frac{1}{f_{\rm
orb}^{1/3}},$ (11)

where $\mathcal{M}_{c}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$ is the system’s
chirp mass, and $D_{L}$ is the systems luminosity distance which we assume to
be the distance of each simulated DWD to the Sun.

For stationary sources, the SNR is modified due to the lack of observable
orbital evolution as

$\rho_{\text{circ,stat}}^{2}=\frac{h_{2}^{2}T_{\text{obs}}}{4S_{n}(f_{2})}$
(12)

with the observation time $T_{\text{obs}}$. Here, $h_{2}$ is strain amplitude
of the source for the second orbital frequency harmonic,

$h_{2}^{2}=\frac{2^{22/3}}{5}\frac{(G\mathcal{M}_{c})^{10/3}}{c^{8}D_{L}^{2}}(\pi
f_{\text{orb}})^{4/3}$ (13)

and is connected to the characteristic strain as

$h_{2}^{2}=\frac{\dot{f}_{2}}{f_{\rm{orb}}^{2}}h_{c}^{2}$ (14)

The amplitude spectral density for a stationary system is finally defined as
$ASD=h_{2}\sqrt{T_{\text{obs}}}$, such that the SNR for stationary source is
simply, $\rho\sim ASD/S_{n}$.

The Galactic foreground included in the Robson et al. (2019) LISA noise curve
was generated using a different binary evolution code and set of model
assumptions for DWD formation and evolution (Toonen et al., 2012; Korol et
al., 2017). Thus, we use the detector curve only and generate an approximate
foreground from each of our populations as follows. Instead of performing a
full source subtraction algorithm (e.g. Littenberg et al., 2020), which is out
of the scope of this work, we calculate the PSD of the Galactic DWD population
with a frequency resolution set by the LISA mission time as
$1/T_{\rm{obs}}\sim 1/4\,\rm{yr}^{-1}\sim 8\times 10^{-9}\,\rm{Hz}$. We then
approximate the foreground as the running median of the PSD with a boxcar
window with a width of $10^{3}$ frequency bins similar to Benacquista &
Holley-Bockelmann (2006). The Galactic DWD PSD is truncated near
$10\,\rm{mHz}$ for both of our models because we remove all DWDs which
experience Roche-lobe overflow. In order to smooth the effect of this
truncation in our foreground, we fit each running median with fourth-order
polynomials for GW frequencies up to $1\,\rm{mHz}$, thus allowing an
approximation of the foreground PSD for higher frequencies. These fits are
listed in Table 1 where the polynomial is described as

$\log_{10}(\rm{confusion\ fit}/\rm{Hz})=a\,x^{4}+b\,x^{3}+c\,x^{2}+d\,x+e$
(15)

and $x=\log_{10}(f_{\rm{GW}}/\rm{Hz})$. We add the fitted polynomial of the
PSD’s running median to the LISA noise PSD to obtain a sensitivity curve for
each model.

model | a | b | c | d | e
---|---|---|---|---|---
F50 | -223.5 | -189.8 | -76.8 | -14.0 | -1.0
FZ | -558.5 | -575.4 | -243.1 | -45.8 | -3.2
Table 1: Polynomial fitting coefficients for the confusion foreground fit of
Equation 15 for each binary fraction model.

##  4 Metallicity effects on the formation and evolution of DWDs

###  4.1 DWD types and their formation channels

As discussed in Section 2.1, we consider four DWD sub-types, which each
contribute differently to LISA’s GW signals: He + He, CO + He, CO + CO, and
ONe + X. Each sub-type has a unique distribution in their formation times,
initial masses, radii, and orbital periods stemming from variations in their
evolution channels and their formation efficiency. Here we describe the
general formation scenarios and population properties of Galactic close DWDs
which may be observable by LISA.

He WDs are unable to form through single star evolution within the lifetime of
the Milky Way. Instead, they originate through interactions in close binary
systems or binaries with large eccentricities. Because of this, He WDs are
able to form with low component masses on order $\sim$ 0.1M⊙, with the
majority of He WDs in our simulations having masses between 0.2-0.5 M⊙. He +
He DWDs form through evolution of close binary systems, during which their
envelopes are both stripped through RLO and CE phase interactions before
Helium ignition occurs. The two progenitor stars generally have masses
$\lesssim 3M_{\odot}$ which is lower than the progenitors of other DWD types.
Our simulated He + He DWDs have an approximately constant distribution of
formation times $\gtrsim$ 2.5 Gyr. Lastly, since the ZAMS separations are
skewed towards shorter values, we also see that the resulting DWD separations
are smaller on average than that of other DWD types.

Figure 2: The DWD formation efficiency vs metallicity of DWD populations
simulated with COSMIC. Each panel shows the formation efficiency for a given
DWD type. The solid lines indicate the formation efficiency for model FZ which
incorporates a metallicity-dependent binary fraction. The dashed lines
indicate the formation efficiency for model F50, which assumes a constant
binary fraction of $50\%$. The DWD formation efficiency drops by a factor of
$4-5$ for model FZ and a factor of 2-5 for model F50. See Section 4.2 for a
careful description of the trends for each DWD type.

A CO WD forms when a star is able to begin the helium burning process before
its envelope is stripped. Thus to form a CO + He DWD, RLO and CE stages occur
after one component experiences core helium burning, but before the other
component can. Most close CO + He DWDs form in approximately 2 Gyr after ZAMS
and with very short periods because the He WD is formed through the ejection
of a common envelope which greatly reduces the orbital separation. Because of
these short formation separations, many CO + He DWDs merge before the present
day. Due to their asymmetric mass distributions, they have lower chirp masses,
but their shorter periods make them important candidates for LISA detection.

To prevent the two stars’ envelopes from being stripped before helium
ignition, CO + CO DWDs typically form from progenitors in wider orbits, and
the two components may have little to no interaction during their evolution
from ZAMS to DWD. CO DWDs thus have a distribution in progenitor separation
that extends to larger values than for other DWD types. Most need $>0.3$ Gyr
to form. These systems have component WD masses between
$0.35$–$1.0\,M_{\odot}$ and make up the majority of the DWD population.

ONe WDs are rare and typically form from massive progenitor stars which evolve
through the asymptotic giant branch (AGB) phase, thus resulting in a higher
mass WD. All ONe WDs in our COSMIC populations, e.g., have progenitor ZAMS
masses above $4\,M_{\odot}$, and the resulting ONe WDs have a relatively flat
distribution of masses from $1.05\,M_{\odot}$ up to the Chandrasekkhar limit
of $1.4\,M_{\odot}$. Because an ONe WD can have a companion of any other WD
type in our study, there is a spread in their distributions for separation,
secondary mass, final orbital period, and formation time. In general, however,
these systems result from wider separations to allow for the evolution of the
ONe component without merging, e.g. all initial separations in our COSMIC
populations have separations $\gtrsim 1.5\,R_{\odot}$. ONe + X DWDs can form
on short timescales, as low as $30\,\rm{Myr}$ for the majority of high-
metallicity systems.

###  4.2 Metallicity-dependent trends in the formation efficiency of DWDs

The number of DWDs formed per unit solar mass of ZAMS star formation, or DWD
formation efficiency $\eta_{\rm{form}}(Z))$, varies with metallicity.
Consequently, a metallicity-dependent binary fraction further impacts the
efficiency of DWD formation within the Galaxy. Figure 2 shows the DWD
formation efficiency as a function of metallicity for each DWD type and for
each of our binary fraction models. In general, the formation efficiency
decreases with increasing metallicity. This effect is exaggerated for model FZ
which assumes a binary fraction which also decreases with increasing
metallicity.

For He + He DWDs, the sudden drop in formation efficiency near
$\log(Z/Z_{\odot})=-1.0$, is generally caused by the timescale for which the
initially more massive star in the DWD progenitor overflows it’s Roche lobe.
At lower metallicities, donors tend to fill their Roche lobes while they are
still on the main sequence and the mass transfer remains stable. This is
because, in our COSMIC models lower metallicity stars evolve faster than high-
metallicity stars and have larger radii near the end of the MS. At higher
metallicities, mass transfer is initiated when the donor has left the main
sequence and the binary enters a common envelope evolution. The ZAMS orbital
period leads to different specific evolutionary channels. For short-period
systems with periods below $10~{}\rm{days}$, this leads to a stellar merger
due to insufficient orbital energy to eject the envelope. Stellar mergers
continue to dominate the evolutionary pathways of systems with intermediate
orbital periods ($1<\log_{10}(P_{\rm{orb}}/\rm{day})<2.5$). However, an
additional growing number of merging systems arises with one WD component and
one stellar companion, and systems which don’t interact at all and thus do not
form a He + He DWD. The distinction between the various scenarios in this
intermediate orbital period range depends on the combination of their ZAMS
masses and orbital periods. Finally, at wider initial periods
$\log_{10}(\rm{P})>2.5$, the decrease in formation efficiency is dominated by
systems which never interact and thus do not form a DWD before the present
day.

Figure 3: Average interaction separation, $\overline{\mathrm{a}}_{\rm{CE}}$ of
progenitors of close DWDs across metallicity for each DWD type from model FZ.
Solid lines show the average separation at the first RLO for binaries in each
metallicity bin. The shaded regions show the $1\sigma$ spread around the mean
within each metallicity bin. The average interaction separation increases with
metallicity for every DWD type. The positive trend in the average interaction
separation is a direct consequence of larger envelope masses of higher-
metallicity donors which are less evolved than their lower-metallicity
counterparts. The sharp drop of the CO + CO average at the highest-metallicity
bin is due to binning effects of the figure: it only contains 115 systems, of
which the majority stem from a single low-separation system which was sampled
94 times times during the simulation.

For short-period CO + He DWD progenitor binaries with orbital periods below
$\sim 30~{}\rm{days}$, the only channel for mergers before DWD formation is
during a CE evolution. This channel is similar to the stellar merger channel
for He + He DWDs, but due to higher-mass progenitors the CE evolution results
in the binary component with a higher mass becoming a CO WD. At higher
metallicities, the CO WD merges with its companion. Higher-metallicity
progenitors have a larger fraction of their mass in the convective envelope
when compared to lower-metallicity stars of the same mass (Amard, L. et al.,
2019; Amard & Matt, 2020) Thus a CE phase with a higher-metallicity progenitor
requires more orbital energy to eject the common envelope, causing a larger
amount of orbital shrinking which results in a merger later in its binary
evolution.

For CO + He DWD progenitor systems with intermediate orbital periods
($1.5\leq\log_{10}(P_{\rm{orb}}/\rm{day})<2.5$) the mechanisms which impact
formation efficiency are complex. The dominant way to lose CO + He DWDs is
stellar mergers that occur during the second CE phase. For lower-metallicity
systems, the second CE is successful, while the opposite is true for higher-
metallicity systems. This is again because the envelope mass of the CE
evolution donor at higher metallicity is larger than it would be at lower
metallicity and thus requires more orbital energy to eject in the CE phase.
This leads to shorter post-CE orbital periods after the first CE phase and
hence mergers during the second phase.

If a lower-metallicity system forms a CO + He DWD and a higher-metallicity
system does not, this is because the initially more massive binary component
initiates a CE while on the giant branch (GB) instead of the asymptotic giant
branch (AGB). This instead leaves behind a He WD with a stellar companion,
which is the dominating scenario that restricts formation efficiency. There
are a few edge cases where either a He + He DWD is formed instead, or when a
CO WD and a stellar companion is formed and there has been very nearly but not
quite enough time for a CO + He DWD to form.

Long-period binaries with $\log_{10}(P_{\rm{orb}}/\rm{day})>2.5$ also display
complex scenarios that hinder CO + He DWD formation. In near-equal
contributions, our COSMIC simulations produce either stellar mergers or stable
non-DWD binaries at the end of the Hubble time. At these orbital periods,
stellar mergers always occur with a CE phase between a stellar companion and a
CO or He WD. Similar to binaries with shorter orbital periods, the mergers
occur because of increased CE donor envelope masses at higher metallicities. A
subdominant channel of stable He + He DWDs can also occur when a high-
metallicity primary overflows it’s Roche lobe while still on the GB and thus
forms a He WD.

The decrease in the CO + CO DWD formation efficiency with increasing
metallicity stems from different evolutionary channels which arise at the ZAMS
orbital period boundary of $\log_{10}(P_{\rm{orb}}/\rm{day})\simeq 3$. We find
that for binaries with orbital periods below this boundary, the most common
way that CO + CO DWDs form at lower metallicities but not at higher
metallicities is through stellar mergers during a CE phase with a donor that
is still on the GB. For lower-metallicity binaries, which evolve on faster
timescales, the primaries enter CE evolution while on the AGB instead and the
binary is able to survive. For binaries with orbital periods above the
boundary, the vast majority of systems with wide initial orbits end up as
stable binaries. The systems which don’t form a CO DWD at high metallicity do
so because one or both of the binary components initiate a CE phase while
still on the GB thus producing a CO + He or He + He DWD.

The strongest effect which hinders formations of higher-metallicity ONe + X
DWDs is the strength of metallicity-dependent stellar winds assumed in our
model (Vink et al., 2001). The strength of line-driven winds varies more
strongly for the more massive ($\geq 5\,\rm{M}_{\odot}$) progenitors of ONe
WDs relative to the other lower mass WD progenitors. At higher metallicities,
ONe WD progenitors can lose enough mass through winds such that they don’t
ignite their CO cores and thus leave behind a CO WD. Conversely, the lower-
metallicity progenitors retain enough mass to cause carbon ignition and leave
behind an ONe WD.

###  4.3 Metallicity trends in DWD progenitor common envelope separation

Figure 4: The number of LISA-band systems formed for each DWD type as a
function of the base-10 logarithm of metallicity, normalized to solar value.
The solid line shows the FZ population with a metallicity-dependent binary
fraction incorporated, and the dashed line shows the F50 population for a
standard binary fraction of 0.5. The LISA population of DWDs is dominated by
stars with super-solar metallicities. This is true even for model FZ, which
drops off significantly for higher metallicities, because of the large number
of stars formed in m12i beyond Z$\simeq$Z⊙. There is a double peak in the He +
He population; the first peak is caused by the sharp drop in formation
efficiency past Z$\simeq$0.1$Z_{\odot}$ which is then greatly overcompensated
for by the amount of star formation at higher metallicities which forms the
second peak.

All systems that end up radiating GWs in the LISA band have undergone at least
one phase of CE evolution. For systems which experience a stable RLO mass
transfer in the first interaction, the CE phase plays a key role in shrinking
systems with initially wide separations to bring them into the LISA band.

Figure 3 shows the average separation at the first instance of CE evolution,
$a_{\text{CE}}$ of all DWD progenitors that result in systems orbiting in the
LISA frequency band at present day, as a function of metallicity for each DWD
type. The solid lines denote the average value, and the $1\sigma$ variance is
shown in the surrounding shading. The average CE separation increases in
general across metallicity. Higher-metallicity binaries will interact earlier
in the binary’s lifetime than a lower-metallicity binary of equal separation
due to their relatively larger maximum radii. Since higher-metallicity
binaries are also more likely to merge during CE interactions because of their
relatively more massive donor envelopes, the DWDs which survive and eventually
orbit in the LISA band originate from systems with higher interaction
separations which allow their orbits to shrink significantly during CE phases
without merging. This regulation plays a key role in smearing out any
observable effects of a metallicity-dependent binary fraction in the
population of DWDs observable by LISA.

Note that the average CE separation is close to the qualitative limit for
close binaries of $P_{\rm{orb}}\gtrsim 10^{4}\,\rm{day}$ for the CO + CO and
ONe + X binaries. We investigated the fraction of close DWD progenitors which
formed with orbital periods above this range for each DWD type. For the He +
He DWDs, less than $10\%$ of all progenitors, regardless of metallicity
originate in wide orbits. For CO + He population, the percentage increases
from $10$–$15\%$ from low to high metallicities, and for the more massive CO +
CO and ONe + X DWDs, the percentages increase from $15$–$25\%$ and $20$–$40\%$
respectively. This suggests that our application of the close binary fraction
may not be valid over the entire range of the DWD progenitor population.
However, since the relative sizes of CO + CO and ONe + X are small compared to
the He + He and CO + He populations, we leave a further investigation of these
effects to a future study.

##  5 Metallicity Dependence of the LISA DWD Population

While metallicity impacts the intrinsic properties of our simulated DWD
populations as described in Sections 4.2 and 4.3, when we consider the
present-day Galactic close DWDs we find that the population detectable by LISA
changes only numerically. The number of DWDs in the LISA frequency band
decreases by $\sim 50\%$ when comparing model F50 to model FZ.

Figure 5: The PSD of the entire Galactic DWD GW foreground, summed over all
metallicities and DWD types vs GW frequency for models FZ and F50. The
vertical lines show the PSD where model F50 is shown in dark blue and model FZ
is shown in light blue. The dashed lines show fits to the rolling boxcar
median of width 1000 bins for each PSD (see 3 for a discussion). A
metallicity-dependent binary fraction (model FZ) yields fewer DWDs across all
frequencies than a $50\%$ binary fraction (model F50) by roughly a factor of
$2$. This produces a lower GW confusion foreground for frequencies with
f${}_{\rm{GW}}\leq\sim 10^{-3}$ Hz.

Figure 6: The ASD vs GW frequency for DWDs resolved with SNR $>7$ for each DWD
type where the top row shows the population from model FZ and the bottom row
shows the population from model F50. In each panel, the LISA sensitivity
curve, including the confusion foreground for each model, is shown in black
and the total population for each model is shown in grey. We find that each
model qualitatively exhibits similar characteristics and that the only change
is in the yield of resolved DWDs for each type based on the strength of the
confusion foreground.

Figure 4 shows the number of DWDs orbiting with frequencies
$f_{\rm{GW}}>0.1\,\rm{mHz}$ against metallicity for each DWD type. The solid
lines show DWDs from model FZ, and the dashed lines show DWDs from model F50.
The He + CO, CO + CO, and ONe + X populations each have strong peaks in the
number of DWDs near solar metallicity at which the majority of star formation
in galaxy m12i occurs. The largest contribution to the population comes from
metallicities above $\sim$ 0.01Z⊙. The discrepancy between the two binary
fraction models is also the most significant above this threshold. When
creating our DWD populations the DWD formation efficiency, number of m12i star
particles, and binary fraction all compete. The amount of stars formed in m12i
at higher metallicity values overwhelms the drop in DWD formation efficiency
by multiple orders of magnitude, so this effect dominates when determining the
number of stars initially sampled in the population for $f_{\rm{b}}$.

There are two peaks in the distribution of He + He DWDs. This occurs because
near $Z\simeq 0.1Z_{\odot}$, the DWD formation efficiency transitions from
near constant values to a sharp decrease (see Figure 2). However, for a drop
in the formation efficiency by a factor of of $\sim$ six, the amount of star
formation in galaxy m12i increases by more than an order of magnitude for
supersolar metallicities. Above $Z=Z_{\odot}$ this overcompensates for the
efficiency drop, producing the second peak.

Figure 7: The chirp mass vs distance for each DWD type is shown. Only systems
with observable evolution in their GW frequency, i.e those which are chirping,
and with SNR $>7$, are shown, since these are systems for which distance can
be separated from chirp mass within their strain amplitude. Each panel shows
one DWD type, summed over all metallicities. Model FZ is indicated with solid
light blue contours, and model F50 is indicated with dark blue dashed contours
respectively. Contours are shown at the $5^{\rm{th}}$, 25th, 50th, 75th, and
95th percentiles. Despite intrinsic changes to population properties induced
by a metallicity-dependent binary fraction, and a reduction in the height of
the DWD Galactic foreground the distributions are very similar.

The reduced number of DWDs in model FZ relative to model F50 is also apparent
in the GW PSD LISA will observe. We show the GW PSD of each model, as well as
the confusion estimate, in Figure 5. While model F50 (dark blue) and model FZ
(light blue) produce several thousand large spikes in the PSD across LISA’s
frequency band, the overall foreground height, including the confusion, is
larger for model F50. This is a direct consequence of the overall reduction in
the size of the close DWD population in model FZ.

Similar to previous studies, we find that LISA will be able to resolve several
thousand DWDs. Figure 6 shows the amplitude spectral density vs GW frequency
of the resolved systems with SNR $>7$ for each DWD type, where the top and
bottom rows show results for models FZ and F50 respectively. For comparison,
we also show the LISA sensitivity curve, including the modeled confusion
foreground from each population in black, and the the entire population for
each model in grey. Apart from each DWD type having a different abundance of
resolved systems, the population-wide characteristics remain unchanged between
the two binary fraction models. The populations containing at least one He WD
occupy the lower-ASD, higher-GW frequency region of parameter space compared
to the total population, with CO + He DWDs having larger ASDs than the He + He
DWD population. Conversely, DWD types without a He WD component tend to occupy
the higher-ASD, lower-GW frequency region of parameter space. This difference
is largely due to the formation scenario of DWDs containing a He WD, which
form from the ejection of a common envelope created by the He WD progenitor.
These lower mass progenitors overflow their Roche lobes at closer separations
relative to the higher mass progenitors (e.g. Figure 3) and thus also produce
closer DWDs. While the distance to any one DWD strongly influences its ASD,
DWD populations without a He WD component have, on average, higher ASDs due to
their more massive WD components.

We note that while the height of the confusion foreground and the number of
DWDs radiating GWs with $f_{\rm{GW}}>10^{-4}\,\rm{Hz}$ is reduced by a factor
of 2 for model FZ relative to model F50, the number of resolved sources is not
reduced to an equal degree. This is because in the absence of less competing
GW signals from DWDs in the LISA frequency band, more DWDs can be individually
resolved. Between models FZ and F50, the number of resolvable systems with SNR
$>7$ over all DWD types decreases by only $\sim 14\%$.

The distance and chirp mass of DWDs which exhibit observable orbital evolution
due to the emission of GWs during the LISA mission can be measured. This is
because the chirp mass – distance degeneracy in the observed strain can be
broken with the observed GW frequency evolution, or chirp. Assuming a chirp
resolution of $1/T_{\rm{obs}}^{2}\sim 8\times 10^{-9}\,\rm{Hz}^{2}$, we select
the DWDs whose chirp masses and distances can be measured. Figure 7 shows the
chirp mass vs the luminosity distance for each DWD type in this selected
population. The contours show the $5^{\rm{th}}$, $25^{\rm{th}}$,
$50^{\rm{th}}$, $75^{\rm{th}}$, and $95^{\rm{th}}$ percentiles for models FZ
(light blue) and F50 (dark blue, dashed). Despite the reduction in the height
of the confusion foreground when considering model FZ relative to F50, we find
that LISA is unable to differentiate between the chirp mass – distance
distributions of the two models.

##  6 Binary evolution assumption variations

Figure 8: The top panel shows the Galactic DWD confusion fit vs GW frequency
for different binary evolution assumption variations (colors) with model F50
shown in dashed lines and model FZ shown in solid lines. The bottom panel
shows the ratio of the number of DWDs orbiting in the LISA frequency band for
model FZ to model F50 vs to number of DWDs orbiting in the LISA band for model
F50 only. While both the height of the confusion foreground and number of LISA
DWDs changes for each variation, the FZ models within each variation exhibit a
constant reduction by a factor of two compared with the F50 models.

In order to test the robustness of the reduction of the height of the Galactic
DWD GW foreground when assuming a metallicity-dependent binary fraction, we
repeat our analysis using three different binary evolution assumption
variations. For each variation, we consider models FZ and F50 as done in the
fiducial case described above. In variation $q3$, we vary the assumption for
the critical mass ratios at which a RLO interaction remains stable or becomes
unstable from our fiducial assumptions such that the critical mass ratio is
increased to $3.0$ and thus allows stable mass transfer for more massive RLO
donors. In variations $\alpha 25$ and $\alpha 5$, we modify the common
envelope ejection efficiency to be either much less ($\alpha=0.25$) or more
($\alpha=5$) than in our fiducial assumption ($\alpha=1$). Larger common
envelope ejection efficiencies lead to wider post-CE separations, while
smaller ejection efficiencies either lead to closer post-CE separations or
mergers where the envelope ejection fails.

In each variation, changing the binary evolution assumptions dramatically
changes the formation and evolution of the DWD populations. These changes lead
to large shifts in the overall number of close DWDs in our synthetic present-
day Milky-Way-like galaxies. We find that the total number of DWDs with
$f_{\rm{GW}}>10^{-4}\,\rm{Hz}$ increases for both variations $q3$ and $\alpha
5$. This is because there are fewer stellar mergers which occur before the
formation of a DWD, thus allowing more systems to evolve due to GW emission
and orbit in the LISA frequency band at present. Conversely, for variation
$\alpha 0.25$, we find that the number of DWDs orbiting with frequencies in
the LISA band is drastically reduced. This is because of the highly
inefficient use of orbital energy to eject the common envelope, thus producing
more stellar mergers, or closer binaries which are more prone to future
mergers.

Interestingly, when we compare the populations of each binary fraction model
for our variations, we find that the number of close DWDs reduces by the
constant factor of two as seen in our fiducial set of assumptions. This is
illustrated in 8. The top panel shows the confusion foreground fits for each
binary evolution variation (different colored lines) and for each binary
fraction model where the solid lines show FZ models and dotted lines show F50
models. The bottom panel shows the ratio of the number of DWDs orbiting in the
LISA frequency band for the FZ models vs the F50 models for each variation.
Even though the number of DWDs in the LISA band spans over two orders of
magnitude, the ratio of FZ to F50, as well as the spectral shape of the
confusion fit, stays fixed at a constant factor of two reduction. This
suggests that assuming a metallicity-dependent binary fraction reduces the
size of the Galactic close DWD population by a factor of $\sim$ two and the
strength of the Galactic DWD GW foreground for LISA _regardless of the chosen
binary evolution model_.

##  7 Conclusions

In this study, we have investigated the effects of assuming a metallicity-
dependent binary fraction on the formation and evolution of the Galactic
population of DWDs with a special focus on the implications for LISA. Based on
our synthetic Milky-Way-like galaxy catalogs of DWDs, we find that applying a
metallicity-dependent binary fraction changes the formation efficiency and
evolutionary history of DWD populations. However, when considering the close
DWD populations observable by LISA, we find that the only distinguishing
features between models which assume a metallicity-dependent binary fraction
(model FZ) and models which assume a flat $50\%$ binary fraction (model F50)
are the population sizes and the strength of the Galactic DWD GW foreground.
Models which assume a metallicity-dependent binary fraction produce Galactic
DWD populations that are reduced by a factor of two relative to the standard
model assumptions. This reduction extends to the height of the confusion
foreground in the LISA data stream.

We extended our study to include three binary evolution assumption variations
to investigate whether the DWD population reduction was robust to changes in
assumptions for mass transfer stability and common envelope ejection
efficiencies. While our binary evolution assumption variations change the size
of the LISA-observable populations dramatically, the reduction in the size of
the close Galactic DWD population and the height of the confusion foreground
for models which assume a metallicity-dependent binary fraction is robust. An
important consequence of a lower Galactic DWD confusion foreground is that
relative to the total DWD population, more DWDs can be individually resolved
because of the reduction in competing GW signals. While the number of DWDs
radiating GWs in the LISA frequency band is reduced by a factor of two for
model FZ relative to model F50, the number of resolved sources is less
affected with a population-wide reduction of $14\%$. These results are far-
reaching since the strength of the Galactic DWD confusion foreground has
direct consequences on the detectability of all other LISA sources with small
SNRs. An increase in resolution capability from the reduced confusion
foreground can be extended to other galactic binaries that LISA will observe
at these frequencies like those involving neutron stars and stellar-origin
black holes, as well as other more exotic GW sources like merging supermassive
black holes, extreme mass ratio inspirals, or cosmological GW backgrounds.
Based on our results, we suggest that studies which employ fits to the
confusion foreground based on population synthesis results consider reducing
the strength of the Galactic foreground PSD by a factor of two.

The authors are grateful for helpful discussions with Carles Badenes,
Christine Mazzola Daher, and the Gravitational Waves and Astronomical Data
groups at the CCA. S.T. was supported by an Undergraduate Student Research
Award (USRA) at CITA from the Natural Sciences and Engineering Research
Council of Canada (NSERC), Reference # 498223. K.B. is grateful for support
from the Jeffrey L. Bishop Fellowship. The Flatiron Institute is supported by
the Simons Foundation.

## Data Availability

All data and software required to reproduce our results are available through
GitHub and Zenodo which are accessible through the icon links associated with
our abstract and each figure.

## References

  * Adams & Cornish (2014) Adams, M. R., & Cornish, N. J. 2014, Phys. Rev. D, 89, 022001, doi: 10.1103/PhysRevD.89.022001
  * Amard & Matt (2020) Amard, L., & Matt, S. P. 2020, The Astrophysical Journal, 889, 108, doi: 10.3847/1538-4357/ab6173
  * Amard, L. et al. (2019) Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77, doi: 10.1051/0004-6361/201935160
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786.  https://arxiv.org/abs/1702.00786
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Babak et al. (2017) Babak, S., Gair, J., Sesana, A., et al. 2017, Phys. Rev. D, 95, 103012, doi: 10.1103/PhysRevD.95.103012
  * Badenes et al. (2018) Badenes, C., Mazzola, C., Thompson, T. A., et al. 2018, ApJ, 854, 147, doi: 10.3847/1538-4357/aaa765
  * Barack & Cutler (2007) Barack, L., & Cutler, C. 2007, Phys. Rev. D, 75, 042003, doi: 10.1103/PhysRevD.75.042003
  * Bartolo et al. (2016) Bartolo, N., Caprini, C., Domcke, V., et al. 2016, J. Cosmology Astropart. Phys, 2016, 026, doi: 10.1088/1475-7516/2016/12/026
  * Bellovary et al. (2019) Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913, doi: 10.1093/mnras/sty2842
  * Benacquista & Holley-Bockelmann (2006) Benacquista, M., & Holley-Bockelmann, K. 2006, ApJ, 645, 589, doi: 10.1086/504024
  * Berti et al. (2006) Berti, E., Cardoso, V., & Will, C. M. 2006, Phys. Rev. D, 73, 064030, doi: 10.1103/PhysRevD.73.064030
  * Boileau et al. (2021) Boileau, G., Lamberts, A., Christensen, N., Cornish, N. J., & Meyer, R. 2021, MNRAS, 508, 803, doi: 10.1093/mnras/stab2575
  * Breivik et al. (2018) Breivik, K., Kremer, K., Bueno, M., et al. 2018, ApJ, 854, L1, doi: 10.3847/2041-8213/aaaa23
  * Breivik et al. (2020a) Breivik, K., Mingarelli, C. M. F., & Larson, S. L. 2020a, ApJ, 901, 4, doi: 10.3847/1538-4357/abab99
  * Breivik et al. (2020b) Breivik, K., Coughlin, S., Zevin, M., et al. 2020b, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Caldwell et al. (2019) Caldwell, R. R., Smith, T. L., & Walker, D. G. E. 2019, Phys. Rev. D, 100, 043513, doi: 10.1103/PhysRevD.100.043513
  * Caprini et al. (2016) Caprini, C., Hindmarsh, M., Huber, S., et al. 2016, J. Cosmology Astropart. Phys, 2016, 001, doi: 10.1088/1475-7516/2016/04/001
  * Chatziioannou et al. (2017) Chatziioannou, K., Clark, J. A., Bauswein, A., et al. 2017, Phys. Rev. D, 96, 124035, doi: 10.1103/PhysRevD.96.124035
  * Claeys et al. (2014) Claeys, J. S. W., Pols, O. R., Izzard, R. G., Vink, J., & Verbunt, F. W. M. 2014, A&A, 563, A83, doi: 10.1051/0004-6361/201322714
  * Cornish (2020) Cornish, N. J. 2020, Phys. Rev. D, 102, 124038, doi: 10.1103/PhysRevD.102.124038
  * de Mink et al. (2007) de Mink, S. E., Pols, O. R., & Hilditch, R. W. 2007, A&A, 467, 1181, doi: 10.1051/0004-6361:20067007
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960
  * Flanagan & Hughes (1998) Flanagan, É. É., & Hughes, S. A. 1998, Phys. Rev. D, 57, 4535, doi: 10.1103/PhysRevD.57.4535
  * Geller et al. (2019) Geller, A. M., Leigh, N. W. C., Giersz, M., Kremer, K., & Rasio, F. A. 2019, ApJ, 872, 165, doi: 10.3847/1538-4357/ab0214
  * Gokhale et al. (2007) Gokhale, V., Peng, X. M., & Frank, J. 2007, ApJ, 655, 1010, doi: 10.1086/510119
  * Goldberg & Mazeh (1994) Goldberg, D., & Mazeh, T. 1994, A&A, 282, 801
  * Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543, doi: 10.1046/j.1365-8711.2000.03426.x
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Jayasinghe et al. (2021) Jayasinghe, T., Kochanek, C. S., Strader, J., et al. 2021, MNRAS, 506, 4083, doi: 10.1093/mnras/stab1920
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kilic et al. (2021) Kilic, M., Bédard, A., & Bergeron, P. 2021, MNRAS, 502, 4972, doi: 10.1093/mnras/stab439
  * Klein et al. (2016) Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003, doi: 10.1103/PhysRevD.93.024003
  * Korol et al. (2021) Korol, V., Hallakoun, N., Toonen, S., & Karnesis, N. 2021, arXiv e-prints, arXiv:2109.10972.  https://arxiv.org/abs/2109.10972
  * Korol et al. (2018) Korol, V., Koop, O., & Rossi, E. M. 2018, ApJ, 866, L20, doi: 10.3847/2041-8213/aae587
  * Korol et al. (2019) Korol, V., Rossi, E. M., & Barausse, E. 2019, MNRAS, 483, 5518, doi: 10.1093/mnras/sty3440
  * Korol et al. (2017) Korol, V., Rossi, E. M., Groot, P. J., et al. 2017, MNRAS, 470, 1894, doi: 10.1093/mnras/stx1285
  * Kremer et al. (2017) Kremer, K., Breivik, K., Larson, S. L., & Kalogera, V. 2017, ApJ, 846, 95, doi: 10.3847/1538-4357/aa8557
  * Kremer et al. (2015) Kremer, K., Sepinsky, J., & Kalogera, V. 2015, ApJ, 806, 76, doi: 10.1088/0004-637X/806/1/76
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Lamberts et al. (2019) Lamberts, A., Blunt, S., Littenberg, T. B., et al. 2019, MNRAS, 490, 5888, doi: 10.1093/mnras/stz2834
  * Leiner & Geller (2021) Leiner, E. M., & Geller, A. 2021, ApJ, 908, 229, doi: 10.3847/1538-4357/abd7e9
  * Littenberg et al. (2020) Littenberg, T. B., Cornish, N. J., Lackeos, K., & Robson, T. 2020, Phys. Rev. D, 101, 123021, doi: 10.1103/PhysRevD.101.123021
  * Mandhai et al. (2021) Mandhai, S., Lamb, G. P., Tanvir, N. R., et al. 2021, arXiv e-prints, arXiv:2109.09714.  https://arxiv.org/abs/2109.09714
  * Marsh et al. (2004) Marsh, T. R., Nelemans, G., & Steeghs, D. 2004, MNRAS, 350, 113, doi: 10.1111/j.1365-2966.2004.07564.x
  * Mazeh et al. (1992) Mazeh, T., Goldberg, D., Duquennoy, A., & Mayor, M. 1992, ApJ, 401, 265, doi: 10.1086/172058
  * Mazzola et al. (2020) Mazzola, C. N., Badenes, C., Moe, M., et al. 2020, MNRAS, 499, 1607, doi: 10.1093/mnras/staa2859
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Moe et al. (2019) Moe, M., Kratter, K. M., & Badenes, C. 2019, ApJ, 875, 61, doi: 10.3847/1538-4357/ab0d88
  * Moore et al. (2017) Moore, C. J., Chua, A. J. K., & Gair, J. R. 2017, Classical and Quantum Gravity, 34, 195009, doi: 10.1088/1361-6382/aa85fa
  * Nelemans et al. (2001) Nelemans, G., Yungelson, L. R., Portegies Zwart, S. F., & Verbunt, F. 2001, A&A, 365, 491, doi: 10.1051/0004-6361:20000147
  * Nissanke et al. (2012) Nissanke, S., Vallisneri, M., Nelemans, G., & Prince, T. A. 2012, ApJ, 758, 131, doi: 10.1088/0004-637X/758/2/131
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, 1.1.1, Zenodo, doi: 10.5281/zenodo.3509134
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224, doi: 10.1103/PhysRev.136.B1224
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Price-Whelan et al. (2020) Price-Whelan, A. M., Hogg, D. W., Rix, H.-W., et al. 2020, ApJ, 895, 2, doi: 10.3847/1538-4357/ab8acc
  * Robson et al. (2019) Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011, doi: 10.1088/1361-6382/ab1101
  * Ruiter et al. (2010) Ruiter, A. J., Belczynski, K., Benacquista, M., Larson, S. L., & Williams, G. 2010, ApJ, 717, 1006, doi: 10.1088/0004-637X/717/2/1006
  * Sanderson et al. (2020) Sanderson, R. E., Wetzel, A., Loebman, S., et al. 2020, ApJS, 246, 6, doi: 10.3847/1538-4365/ab5b9d
  * Sepinsky & Kalogera (2014) Sepinsky, J. F., & Kalogera, V. 2014, ApJ, 785, 157, doi: 10.1088/0004-637X/785/2/157
  * Shen (2015) Shen, K. J. 2015, ApJ, 805, L6, doi: 10.1088/2041-8205/805/1/L6
  * Soberman et al. (1997) Soberman, G. E., Phinney, E. S., & van den Heuvel, E. P. J. 1997, A&A, 327, 620.  https://arxiv.org/abs/astro-ph/9703016
  * Stephan et al. (2019) Stephan, A. P., Naoz, S., Ghez, A. M., et al. 2019, ApJ, 878, 58, doi: 10.3847/1538-4357/ab1e4d
  * Toonen et al. (2012) Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70, doi: 10.1051/0004-6361/201218966
  * Tout et al. (1997) Tout, C. A., Aarseth, S. J., Pols, O. R., & Eggleton, P. P. 1997, MNRAS, 291, 732, doi: 10.1093/mnras/291.4.732
  * Valsecchi et al. (2012) Valsecchi, F., Farr, W. M., Willems, B., Deloye, C. J., & Kalogera, V. 2012, ApJ, 745, 137, doi: 10.1088/0004-637X/745/2/137
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574, doi: 10.1051/0004-6361:20010127
  * Wagg et al. (2021) Wagg, T., Breivik, K., & de Mink, S. E. 2021, arXiv e-prints, arXiv:2111.08717.  https://arxiv.org/abs/2111.08717
  * Wang et al. (2021) Wang, H., Stephan, A. P., Naoz, S., Hoang, B.-M., & Breivik, K. 2021, ApJ, 917, 76, doi: 10.3847/1538-4357/ac088d
  * Webbink (1985) Webbink, R. F. 1985, Stellar evolution and binaries, ed. J. E. Pringle & R. A. Wade, 39
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Wetzel et al. (2016) Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJ, 827, L23, doi: 10.3847/2041-8205/827/2/L23
  * Wong et al. (2021) Wong, K. W. K., Breivik, K., Kremer, K., & Callister, T. 2021, Phys. Rev. D, 103, 083021, doi: 10.1103/PhysRevD.103.083021
  * Yu & Jeffery (2013) Yu, S., & Jeffery, C. S. 2013, MNRAS, 429, 1602, doi: 10.1093/mnras/sts445
  * Zevin et al. (2020a) Zevin, M., Kelley, L. Z., Nugent, A., et al. 2020a, ApJ, 904, 190, doi: 10.3847/1538-4357/abc266
  * Zevin et al. (2020b) Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020b, ApJ, 899, L1, doi: 10.3847/2041-8213/aba74e
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152, doi: 10.3847/1538-4357/abe40e