# Gravitational wave sources in our Galactic backyard:
Predictions for BHBH, BHNS and NSNS binaries detectable with LISA

T. Wagg Department of Astronomy, University of Washington, Seattle, WA, 98195
Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA  Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße
1, 85741 Garching, Germany  F.S. Broekgaarden Center for Astrophysics —
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  S.E. de
Mink Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany  Anton Pannekoek Institute for Astronomy & Grappa,
University of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands
Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA  N. Frankel Canadian Institute for Theoretical Astrophysics,
University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada  Max
Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany
David A. Dunlap Department for Astronomy and Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada  L.A.C. van Son
Center for Astrophysics — Harvard & Smithsonian, 60 Garden Street, Cambridge,
MA 02138, USA  Anton Pannekoek Institute for Astronomy & Grappa, University of
Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands  Max-Planck-
Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany
S. Justham School of Astronomy & Space Science, University of the Chinese
Academy of Sciences, Beijing 100012, China  Anton Pannekoek Institute for
Astronomy & Grappa, University of Amsterdam, Postbus 94249, 1090 GE Amsterdam,
The Netherlands  Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-
Straße 1, 85741 Garching, Germany Tom Wagg<EMAIL_ADDRESS>
###### Abstract

Future searches for gravitational waves from space will be sensitive to double
compact objects (DCOs) in our Milky-Way. We present new simulations of the
populations of double black holes (BHBHs), black hole neutron stars (BHNSs)
and double neutron stars (NSNSs) that will be detectable by the planned space-
based gravitational wave detector LISA. For our estimates, we use an
empirically-informed model of the metallicity dependent star formation history
of the Milky Way. We populate using an extensive suite of binary population-
synthesis predictions for varying assumptions relating to mass transfer,
common-envelope, supernova kicks, remnant masses and wind mass loss physics.

For a 4(10)-year LISA mission, we predict between 30-370(50-550) detections
over these variations, out of which 6-154(9-238) are BHBHs, 2-198(3-289) are
BHNSs and 3-35(4-57) are NSNSs. We discuss how the variations in the physics
assumptions alter the distribution of properties of the detectable systems,
even when the detection rates are unchanged. In particular we discuss the
observable characteristics such as the chirp mass, eccentricity and sky
localisation and how the BHBH, BHNS and NSNS populations can be distinguished,
both from each other and from the more numerous double white dwarf population.
We further discuss the possibility of multi-messenger observations of pulsar
populations with the Square Kilometre Array (SKA) and assess the benefits of
extending the LISA mission.

gravitational waves, gravitational wave detectors, compact objects, stellar
mass black holes, neutron stars, binary stars, stellar evolution, pulsars

††software: We used LEGWORK to evolve sources over time and calculate signal-
to-noise ratios (Wagg et al., 2021). It is freely available at
https://legwork.readthedocs.io/en/latest/. Simulations in this paper made use
of the COMPAS rapid binary population synthesis code (Team COMPAS: J. Riley et
al., 2021). The simulations performed in this work were simulated with a
COMPAS version that predates the publicly available code. Our version is most
similar to v02.13.01 of the publicly available COMPAS code. Requests for the
original code can be made to Floor Broekgaarden. The authors used STROOPWAFEL
from (Broekgaarden et al., 2019), publicly available at
https://github.com/FloorBroekgaarden/STROOPWAFEL666For the latest pip
installable version of STROOPWAFEL please contact Floor Broekgaarden.. The
authors also made use of Python (v3.8), available at http://www.python.org. In
addition the following Python packages were used: matplotlib (Hunter, 2007),
NumPy (Harris et al., 2020), Astropy (Astropy Collaboration et al., 2013,
2018), Seaborn (Waskom, 2021), SciPy (Virtanen et al., 2020), h5py (Collette
et al., 2021) and Jupyter Lab (Kluyver et al., 2016). This research has made
use of NASA’s Astrophysics Data System Bibliographic Services. We also made
use of the computational facilities from the Harvard FAS Research Computing
cluster.

##  1 Introduction

Since the first direct observation of gravitational waves (Abbott et al.,
2016), the number of black hole (BH) and neutron star (NS) binaries observed
by ground-based gravitational-wave detectors has rapidly grown (Abbott et al.,
2019, 2020b; The LIGO Scientific Collaboration et al., 2021c), offering
exciting insights into the formation, lives and deaths of massive (binary)
stars (e.g. Abbott et al., 2021).

The Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al., 2017; Colpi
et al., 2019) will provide observations in an entirely new regime of
gravitational waves. LISA will observe at lower frequencies ($10^{-5}\lesssim
f/\,\mathrm{Hz}\lesssim 10^{-1}$) than ground-based detectors and so will
enable the study of sources that are imperceptible by ground-based detectors,
such as the mergers of supermassive black holes and extreme mass-ratio
inspirals (e.g. Begelman et al., 1980; Klein et al., 2016). Moreover, this
frequency regime is also of interest for the detection of local stellar-mass
double compact objects (DCOs) millions of years before their merger. This
presents an opportunity for both multi-messenger detections to search for
electromagnetic counterparts, as well as multiband gravitational-wave
detections that can help to constrain binary characteristics (e.g. Sesana,
2016; Gerosa et al., 2019). In addition, LISA will be able to measure the
eccentricities of DCOs, which may yield further constraints on binary
evolution, differentiate between formation channels and distinguish between
DCO types (e.g. Nelemans et al., 2001; Breivik et al., 2016; Antonini et al.,
2017; Rodriguez et al., 2018). Unlike ground-based detectors, LISA only
detects stellar-mass sources in local galaxies, with the majority residing in
the Milky Way. These sources could be used as a probe for the structure of our
galaxy (e.g. Korol et al., 2019).

Traditionally, predictions about the detection of stellar-mass sources with
LISA focus on double white dwarf (WDWD) binaries, as they are abundantly
present in our galaxy and are expected to be the dominant source of stellar-
mass binaries that are detectable by LISA (Nelemans et al., 2001; Ruiter et
al., 2010; Yu & Jeffery, 2010; Nissanke et al., 2012; Korol et al., 2017;
Lamberts et al., 2018). More recently, interest has grown in the detection of
NS and BH binaries. Although they are more rare, LISA detections of these
sources are potentially valuable for learning more about the evolution and
endpoints of massive stars. In this paper we focus on making LISA predictions
for double black hole binaries (BHBHs), black hole neutron star binaries
(BHNSs) and double neutron star binaries (NSNSs).

The detection of NSNSs in LISA could improve our understanding of many
phenomena. Galactic NSNSs have been observed with electromagnetic signals for
several decades (e.g. Hulse & Taylor, 1975; Antoniadis et al., 2016, see also
refs. in Tauris et al., 2017) and more recently the mergers of NSNSs have been
detected with ground-based gravitational-wave detectors (e.g. Abbott et al.,
2017a). A LISA detectable NSNS with a pulsar component close to merger would
be ideal for connecting these populations, as the binary could be observed
from inspiral to merger. NSNSs (and possibly BHNSs) are useful sources for
understanding the origin of r-process elements (e.g. Eichler et al., 1989) as
well as the electromagnetic counterparts to gravitational-wave signals, such
as kilonovae (e.g. Li & Paczyński, 1998; Metzger, 2017), short gamma-ray
bursts (e.g. Berger, 2014), radio emission (e.g. Hotokezaka et al., 2016) and
neutrinos (e.g. Kyutoku et al., 2018).

BHBHs in the Milky Way present a greater observational challenge. To date, no
BH has been observed in a BHBH binary in the Milky Way, and so LISA could
provide the first detection of a Galactic BHBH. The only confirmed BHs in our
galaxy have been discovered as components of X-ray binaries with companion
stars (e.g. Bolton, 1972; Webster & Murdin, 1972). These detections have
observed BHs with masses mainly constrained between $5$ and
$10\,\mathrm{M_{\odot}}$ (Corral-Santana et al., 2016), a stark contrast to
the more massive BHs observed with LIGO/Virgo that tend to contain at least
one BH with a mass greater than $10\,\mathrm{M_{\odot}}$ (Abbott et al.,
2020b). These observations of X-ray binaries suggest the presence of a lower
mass gap (from $2$-$5\,\mathrm{M_{\odot}}$) in which there are no strong
candidates for either black holes or neutron stars (Özel et al., 2010; Farr et
al., 2011) but the gap’s existence remains an open question (e.g. Kreidberg et
al., 2012; Mandel & Müller, 2020). Recently there has also been increased
discussion over the maximum BH mass in our galaxy, with the claims of a
$70\,\mathrm{M_{\odot}}$ BH (Liu et al., 2019) which has subsequently been
challenged (El-Badry & Quataert, 2020; Abdul-Masih et al., 2020; Shenar et
al., 2020; Eldridge et al., 2020, see also Liu et al., 2020) and revised
measurements of the mass of Cygnus X-1 (Miller-Jones et al., 2021). A sample
of BHBHs detected with LISA could possibly help to constrain the BH mass
distribution.

One particularly interesting class of potential LISA sources is BHNSs. With
the recent detection of two BHNSs by the LIGO scientific collaboration, the
existence of these DCOs has been confirmed (The LIGO Scientific Collaboration
et al., 2021b). However, with only two detections (not including the low-
confidence candidate GW190426, Abbott et al. 2020b, or GW190425, GW190814 and
GW190917 which have not been ruled out as BHNSs, Abbott et al. 2020a, c; The
LIGO Scientific Collaboration et al. 2021a) and no electromagnetic
counterparts, the formation rate and properties of BHNSs are still uncertain.
Current predictions for the merger rate of BHNSs range across three orders of
magnitude (e.g. Abadie et al., 2010; Mandel & Broekgaarden, 2021) so the
number of detections in LISA will be important in reducing this uncertainty,
thereby refining our understanding of the remnants and evolution of massive
stars. Similar to NSNSs, these binaries are also expected to have
electromagnetic counterparts. A distinctly exciting possibility is the
detection of a pulsar–BH system or millisecond pulsar–BH system (Narayan et
al., 1991; Pol et al., 2021). These systems could be observed not only by
LISA, but also radio telescopes such as MeerKAT and the Square Kilometre Array
(SKA, Dewdney et al. 2009), which would help to improve the measurement of
individual system parameters and to constrain uncertain binary evolution
processes (e.g. Pfahl et al., 2005; Chattopadhyay et al., 2020).

For the purposes of this investigation, we consider the ‘classical’ isolated
binary evolution channel (e.g. Tutukov & Yungelson, 1973, 1993; Smarr &
Blandford, 1976; Srinivasan, 1989; Kalogera et al., 2007; Belczynski et al.,
2016) in which double compact objects are formed following common-envelope
ejection or a phase of highly non-conservative mass transfer (van den Heuvel,
2011; van den Heuvel et al., 2017). We do not, however, account for several
alternative proposed formation channels, which could affect the rate and
distribution of detectable NS and BH binaries in LISA. These channels include:
dynamical formation in dense star clusters (e.g. Sigurdsson & Hernquist, 1993;
Portegies Zwart & McMillan, 2000; Miller & Lauburg, 2009; Rodriguez et al.,
2015), young/open star clusters (e.g. Ziosi et al., 2014; Di Carlo et al.,
2020; Banerjee, 2020; Rastello et al., 2020, 2021) and (active) galactic
nuclei discs (e.g. Morris, 1993; Antonini & Rasio, 2016; McKernan et al.,
2020), isolated (hierarchical) triple evolution involving Kozai-Lidov
oscillations (e.g. Stephan et al., 2016; Silsbee & Tremaine, 2017; Antonini et
al., 2017; Toonen et al., 2020), and chemically homogenous evolution through
efficient rotational mixing (e.g. de Mink et al., 2009; Mandel & de Mink,
2016; Marchant et al., 2016, 2017; du Buisson et al., 2020).

In this paper, we present models for the detection rate and distribution of
binary properties (masses, frequency, eccentricity, distance, merger time) of
BHBHs, BHNSs and NSNSs formed through isolated binary evolution in the Milky
Way. We explore the effect of varying physical assumptions in our population
synthesis model on our results as well as discuss the effect of extending the
LISA mission length and the prospects for distinguishing DCO detections from
the WDWD background.

Earlier work on BHBHs, BHNSs and NSNSs in LISA has used a variety of
population synthesis codes, Milky Way models and LISA specifications,
resulting in a wide range of predictions (Nelemans et al., 2001; Belczynski et
al., 2010; Liu & Zhang, 2014; Lamberts et al., 2019; Lau et al., 2020; Breivik
et al., 2020; Sesana et al., 2020; Shao & Li, 2021). We build upon previous
efforts but with several important improvements. We explore the effects of
varying binary physics assumptions by repeating our analysis for 20 different
models and comparing the effects on the detection rate and distributions of
source parameters (Broekgaarden et al. 2021, Broekgaarden et al. (in prep.)).
We use a model for the Milky Way that accounts for the chemical enrichment
history and is calibrated on the APOGEE survey (Majewski et al., 2017; Frankel
et al., 2018), whereas most others did not consider the effect of metallicity
in detail (see however Lamberts et al., 2019; Sesana et al., 2020). We provide
a detailed treatment of the eccentricity of detectable sources, both for the
inspiral evolution as well as gravitational wave signal during the LISA
mission. Moreover, the grid of binary population synthesis simulations that we
use is the most extensive of its kind to date and makes use of the adaptive
sampling algorithm STROOPWAFEL (Broekgaarden et al., 2019, 2021). Overall over
2 billion massive binaries were simulated to produce the DCO populations used
in this work. We find that this large number of simulations is important to
reduce the sampling noise even when using adaptive importance sampling.

All data related to the predictions made in this study are publicly available
on Zenodo at Wagg et al. (2021), as are the populations used in our
simulations at Broekgaarden (2021a, BHBH) Broekgaarden (2021b, BHNS) and
Broekgaarden (2021c, NSNS). We make all code used to produce our results
available in a Github repository 111https://github.com/TomWagg/detecting-DCOs-
in-LISA. In addition, the repository contains step-by-step Jupyter notebooks
that explain how to reproduce and change each figure in the paper. In a
companion paper, Wagg et al. (2021), we present
LEGWORK222https://legwork.readthedocs.io, the LISA Evolution and Gravitational
Wave Orbit Kit, a python package designed for making predictions for the
detection of sources with LISA, which we use in this work.

Our paper is structured as follows. In Section 2, we describe the methods for
synthesising a population of binaries, the variations of physical assumptions
that we consider, how we simulate the Milky Way distribution of DCOs and our
methods for calculating a detection rate for LISA. We present the main results
for our fiducial model in Section 3, before exploring the variations in the
detectable population when changing our physical assumptions in Section 4. In
Section 5 we discuss these results. In Section 6, we compare and contrast our
methods and findings to previous work and finish with our conclusions in
Section 7.

##  2 Method

To produce predictions for the DCOs that are detectable with LISA, we use a
synthesised population of DCOs, simulated using the methods described in
Section 2.1. In Section 2.2 we describe our model for the Milky Way and how we
place DCOs in randomly sampled Milky Way instances. We evolve the orbit of
each DCO in a Milky Way instance up to the LISA mission and calculate the
detection rate for that instance using the methods presented in Section 2.3.

###  2.1 Binary population synthesis

We use the grid of 20 binary population synthesis simulations recently
presented in Broekgaarden et al. (2021) and Broekgaarden et al. (in prep.).
This grid of simulations is synthesised using the rapid population synthesis
code COMPAS (Stevenson et al., 2017; Vigna-Gómez et al., 2018; Stevenson et
al., 2019; Broekgaarden et al., 2019). COMPAS follows the approach of the
population synthesis code BSE (Hurley et al., 2000, 2002) and uses fitting
formula and rapid algorithms to efficiently predict the final fate of binary
systems. The code is open source and documented in the papers listed above,
the online documentation333https://compas.science and in the methods paper
(Team COMPAS: J. Riley et al., 2021). We summarise the main assumptions and
settings relevant for this work in Appendix A.

The result of the simulations is a sample of binaries, which, for each
metallicity $Z$, have $N_{\rm binaries}$ binaries with parameters

$\mathbf{b}_{{Z,i}}=\\{m_{1},m_{2},a_{\rm DCO},e_{\rm DCO},t_{\rm
evolve},t_{\rm inspiral},w\\},$ (1)

for $i=1,2,\dots,N_{\rm binaries}$, where $m_{1}$ and $m_{2}$ are the primary
and secondary masses, $a_{\rm DCO}$ and $e_{\rm DCO}$ are the semi-major axis
and eccentricity at the moment of double compact object (DCO) formation,
$t_{\rm evolve}$ is the time between the binary’s zero-age main sequence and
DCO formation, $t_{\rm inspiral}$ is the time between DCO formation (that is,
immediately after the second supernova in the system) and gravitational-wave
merger and $w$ is the adaptive importance sampling weight assigned by
STROOPWAFEL (Broekgaarden et al., 2019, Eq. 7). We sample from these sets of
parameters when creating synthetic galaxies.

###  2.2 Galaxy synthesis

In order to estimate a detection rate of DCOs with statistical uncertainties,
we create a series of random instances of the Milky Way, each populated with a
subsample drawn (with replacement) from the synthesised binaries described in
Section 2.1.

Most previous studies that predict a detection rate for LISA place binaries in
the Milky Way independently of their age or evolution. We improve upon this as
the first study to use an empirically-informed analytical model of the Milky
Way that takes into account the galaxy’s enrichment history by applying the
metallicity-radius-time relation from Frankel et al. (2018). Those authors
developed this relation in order to measure the global efficiency of radial
migration in the Milky Way and calibrated it using a sample of red clump stars
measured with APOGEE (Majewski et al., 2017).

In Section 2.2.1, we outline our model for the Milky Way and in Section 2.2.2
we explain how we combine our population of synthesised DCOs with this Milky
Way model.

####  2.2.1 Milky Way model

Figure 1: A schematic illustrating how we model the Milky Way. The left panel
illustrates the different model aspects: star formation history of three
galactic components (individually shown in the dotted lines), radial
distribution, metallicity-radius-time relation, and height distribution. The
right panel shows an example instance of the Milky Way with $250000$ binaries
shown as points coloured by metallicity. The top panel shows a side-on view
and the bottom panel a face-on view.  .

Fig. 1 shows the distributions and relations outlined in this section and also
displays an example random galaxy drawn using this model.

Our model for the Milky Way accounts for the
low-$[\alpha/\mathrm{Fe}]$444Nomenclature used to describe the enhancement of
$\alpha$ elements compared to iron in stellar atmospheres disc,
high-$[\alpha/\mathrm{Fe}]$ disc and a central component approximating a
bar/bulge. The low- and high-$[\alpha/\mathrm{Fe}]$ discs are often also
referred to as the thin and thick discs because the stellar vertical
distribution is better fit by a double exponential rather than a single one.
However, this doesn’t allow one to assign a star to either the thin or thick
disk purely based on its height above the Galactic plane. Therefore, we
instead use the chemical definition of the two disks (applying the
$[\alpha/\mathrm{Fe}]$ nomenclature) as there is a clear bimodal distribution
in the chemical plane, allowing stars to be assigned to each of the disc
components based on their chemical abundances. For each of the three
components, we use a separate star formation history and spatial distribution,
which we combine into a single model, weighting each component by its present-
day stellar mass. Licquia & Newman (2015) gives that the stellar mass of the
bulge is $0.9\times 10^{10}\,\mathrm{M_{\odot}}$ and the stellar mass of the
disc is $5.2\times 10^{10}\,\mathrm{M_{\odot}}$, which we split equally
between the low- and high-$[\alpha/\mathrm{Fe}]$ discs (e.g., Snaith et al.,
2014).

Star formation history: We use an exponentially declining star formation
history (Frankel et al., 2018) (inspired by the average cosmic star formation
history) for the combined low- and high-$[\alpha/\mathrm{Fe}]$ discs,

$p(\tau)\propto\exp\quantity(-\frac{(\tau_{m}-\tau)}{\tau_{\rm SFR}}),$ (2)

where $\tau$ is the lookback time (the amount of time elapsed between the
binary’s zero-age main sequence and today), $\tau_{m}=12\,\mathrm{Gyr}$ is the
assumed age of the Milky Way and $\tau_{\rm SFR}=6.8\,\mathrm{Gyr}$ is the
star formation timescale (Frankel et al., 2018). The two discs form stars in
mutually exclusive time periods, such that the high-$[\alpha/\mathrm{Fe}]$
disc forms stars in the early history of the galaxy ($8$–$12\,\mathrm{Gyr}$
ago) and the low-$[\alpha/\mathrm{Fe}]$ disc forms stars more recently
($0$–$8\,\mathrm{Gyr}$ ago). Both distributions are normalised so that an
equal amount of mass is formed in each of the two components over their
respective star forming periods.

The star formation history of the Milky Way bulge (which we assume here to be
dominated by the central bar) has many uncertainties due to the (1) sizeable
age measurement uncertainties at large ages in observational studies, (2)
complex selection processes affecting the observed age distributions, and (3)
formation mechanisms that are still under debate. However, the central bar was
shown to contain stars with an extended age range, with most observed stars
between $6$ and $12\,\mathrm{Gyr}$ with a younger tail of ages that could come
from the subsequent secular growth of the Galactic bar (e.g., Bovy et al.,
2019). To model the bar’s age distribution more realistically than in previous
studies (which assume an old bulge coming from a single starburst), we choose
to adopt a more extended star formation history using a $\beta(2,3)$
distribution, shifted and scaled such that stars are only formed in the range
$[6,12]\,\mathrm{Gyr}$. We show these distributions in the top left panel of
Fig. 1.

Radial distribution: For each of the three components we employ the same
single exponential distribution (but with different scale lengths)

$p(R)=\exp(-\frac{R}{R_{d}})\frac{R}{R_{d}^{2}},$ (3)

where $R$ is the Galactocentric radius and $R_{d}$ is the scale length of the
component. For the low-$[\alpha/\mathrm{Fe}]$ disc, we set $R_{d}=R_{\rm
exp}(\tau)$, where $R_{\rm exp}(\tau)$ is the scale length presented in
Frankel et al. (2018, Eq. 6)

$R_{\rm exp}(\tau)=4\,\mathrm{kpc}\quantity(1-\alpha_{R_{\rm
exp}}\quantity(\frac{\tau}{8\,\mathrm{Gyr}})),$ (4)

where $\alpha_{R_{\rm exp}}=0.3$ is the inside-out growth parameter555We find
that $R_{\rm exp}(\tau)=4$ kpc fits the data well and adopt this value rather
than the 3 kpc quoted in Frankel et al. (2018), which was a fixed parameter
(not a fit).

This scale length accounts for the inside-out growth of the
low-$[\alpha/\mathrm{Fe}]$ disc and hence is age dependent. We assume
$R_{d}=(1/0.43)\,\mathrm{kpc}$ for the high-$[\alpha/\mathrm{Fe}]$ disc (Bovy
et al., 2016, Table 1) and $R_{d}=1.5\,\mathrm{kpc}$ for the bar component
(Bovy et al., 2019). We show the combination of these distributions in the
second panel on the left in Fig. 1.

Vertical distribution: Similar to the radial distribution, we use the same
single exponential distribution (but with different scale heights) for each
component, given by

$p(\absolutevalue{z})=\frac{1}{z_{d}}\exp\quantity(-\frac{z}{z_{d}}),$ (5)

where $z$ is the vertical displacement above the Galactic plane and $z_{d}$ is
the scale height. We set $z_{d}=0.3\,\mathrm{kpc}$ for the
low-$[\alpha/\mathrm{Fe}]$ disc (McMillan, 2011) and
$z_{d}=0.95\,\mathrm{kpc}$ for the high-$[\alpha/\mathrm{Fe}]$ disc (Bovy et
al., 2016). For the bar, we set $z_{d}=0.2\,\mathrm{kpc}$ (Wegg et al., 2015).
We show the combination of these distributions in the bottom left panel of
Fig. 1.

Metallicity-radius-time relation: To account for the chemical enrichment of
star forming gas as the Milky Way evolves, we adopt the relation given by
(Frankel et al., 2018, Eq. 7)

$\begin{split}[{\rm Fe}/{\rm H}](R,\tau)&=F_{m}+\nabla[{\rm Fe}/{\rm H}]R\\\
&-\quantity(F_{m}+\nabla[{\rm Fe}/{\rm H}]R^{\rm now}_{[{\rm Fe}/{\rm
H}]=0})f(\tau),\end{split}$ (6)

where

$f(\tau)=\quantity(1-\frac{\tau}{\tau_{m}})^{\gamma_{[{\rm Fe}/{\rm H}]}},$
(7)

$F_{m}=-1\,\mathrm{dex}$ is the metallicity of the gas at the center of the
disc at $\tau=\tau_{m}$, $\nabla[{\rm Fe}/{\rm H}]=-0.075\,\mathrm{kpc^{-1}}$
is the metallicity gradient, $R^{\rm now}_{[{\rm Fe}/{\rm
H}]=0}=8.7\,\mathrm{kpc}$ is the radius at which the present day metallicity
is solar and $\gamma_{[{\rm Fe}/{\rm H}]}=0.3$ sets the time dependence of the
chemical enrichment. We convert this to the representation of metallicity that
we use in this paper by applying (e.g. Bertelli et al., 1994)

$\log_{10}(Z)=0.977[{\rm Fe}/{\rm H}]+\log_{10}(Z_{\odot}).$ (8)

Although Frankel et al. (2018) only fit this model for the
low-$[\alpha/\mathrm{Fe}]$ disc, we also use this metallicity-radius-time
relation for the high-$\alpha$ disc and the bar, but focusing on the chemical
tracks more representative to the inner disc and large ages. Sharma et al.
(2020) showed that using a simple continuous model for both the low- and
high-$[\alpha/\mathrm{Fe}]$ discs, the Milky Way abundance distributions could
be well reproduced. Empirically, the abundance tracks in the
[$\alpha$/Fe]-[Fe/H] plane (and other elements) of the stars in the bulge/bar
follow the same track as those of the old stars in the Solar neighbourhood
(Griffith et al., 2021; Bovy et al., 2019, Fig. 7,), which motivates our
modelling choice to use the same metallicity-radius-time relation.

####  2.2.2 Combining population and galaxy synthesis

For each Milky Way instance, we randomly sample the following set of
parameters

$\mathbf{g}_{{j}}=\\{\tau,R,Z,z,\theta\\}$ (9)

for $j=1,2,\dots,N_{\rm MW}$, where we set $N_{\rm MW}=2\times 10^{5}$,
$\tau,R,Z$ and $z$ are defined and sampled using the distribution functions
specified in Section 2.2.1, $\theta$ is the azimuthal angle sampled uniformly
on $[0,2\pi)$ and $Z$ is the metallicity. Fig. 1 shows an example of a random
Milky Way instance created with these distributions. This shows how these
distributions translate to positions and illustrates the gradient in
metallicity over radius.

We match each set of galaxy parameters $\mathbf{g}_{{j}}$, to a random set of
binary parameters $\mathbf{b}_{{Z,i}}$, by drawing a binary from the closest
metallicity bin to the metallicity in $\mathbf{g}_{{j}}$. If the metallicity
in $\mathbf{g}_{{j}}$ is below the minimum COMPAS metallicity bin
($Z=10^{-4}$), we use this minimum bin. If the metallicity in
$\mathbf{g}_{{j}}$ is above the maximum COMPAS metallicity bin ($Z=0.03$), we
use a randomly selected bin from the five highest metallicity bins.

Each binary is likely to move from its birth orbit. Although all stars in the
Galactic disc experience radial migration (Sellwood & Binney, 2002; Frankel et
al., 2018), DCOs generally experience stronger dynamical evolution as a result
of the effects of both Blaauw kicks (Blaauw, 1961) and natal kicks (e.g. Hobbs
et al., 2005).

The magnitude of the systemic kicks are typically small compared to the
initial circular velocity of a binary at each Galactocentric radius.
Therefore, we expect that kicks will not significantly alter the overall
distribution of their positions (see however, e.g., Brandt & Podsiadlowski
1995; Abbott et al. 2017b). Given this, and for the sake of computational
efficiency, we do not account for the displacement due to systemic kicks in
our analysis.

###  2.3 Gravitational wave detection

We use the Python package LEGWORK (Wagg et al., 2021) to evolve binaries and
calculate their LISA detectability. For a full derivation of the equations
given below see (Wagg et al., 2021, Section 3), or the LEGWORK documentation .

####  2.3.1 Inspiral evolution

Each binary loses orbital energy to gravitational waves throughout its
lifetime. This causes the binary to shrink and circularise over time. In order
to assess the detectability of a binary, we need to know its eccentricity and
frequency at the time of the LISA mission. For each binary in our simulated
Milky Way, we know that the time from DCO formation to today is $\tau-t_{\rm
evolve}$ and that the initial eccentricity and semi-major axis are $e_{\rm
DCO}$ and $a_{\rm DCO}$. We find the eccentricity of the binary at the start
of the LISA mission, $e_{\rm LISA}$, by numerically integrating its time
derivative (Peters, 1964, Eq. 5.13) given the initial conditions. This can be
converted to the semi-major axis at the start of LISA, $a_{\rm LISA}$(Peters,
1964, Eq. 5.11), which in turn gives the orbital frequency, $f_{\rm
orb,LISA}$, by Kepler’s third law since we know the component masses.

####  2.3.2 Binary detectability

We define a binary as detectable if its gravitational wave signal has a
signal-to-noise ratio (SNR) of greater than 7 (e.g. Breivik et al., 2020;
Korol et al., 2020). The sky-, polarisation- and orientation-averaged signal-
to-noise ratio, $\rho$, of an inspiraling binary can be calculated with the
following (e.g. Finn & Thorne, 2000)

$\rho^{2}=\sum_{n=1}^{\infty}\int_{f_{n,i}}^{f_{n,f}}\frac{h_{c,n}^{2}}{f_{n}^{2}S_{\rm
n}\left(f_{n}\right)}\differential{f_{n}},$ (10)

where $n$ is a harmonic of the gravitational wave signal, $f_{n}=n\cdot f_{\rm
orb}$ is the frequency of the $n^{\rm th}$ harmonic of the gravitational wave
signal, $f_{\rm orb}$ is the orbital frequency, $S_{\rm n}(f_{n})$ is the LISA
sensitivity curve at frequency $f_{n}$ (e.g. Robson et al., 2019) and
$h_{c,n}$ is the characteristic strain of the $n^{\rm th}$ harmonic, given by
(e.g. Barack & Cutler, 2004)

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

where $D_{L}$ is the luminosity distance to the source, $f_{\rm orb}$ is the
orbital frequency, $g(n,e)$ and $F(e)$ are given in Peters & Mathews (1963)
and $\mathcal{M}_{c}$ is the chirp mass, defined as

$\mathcal{M}_{c}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}.$ (12)

Note that increasing the length of the LISA mission allows more time for a DCO
to evolve over the mission. Therefore the frequency limits in Eq. 10 are
dictated by the LISA mission length. The SNR generally scales as $\sqrt{T_{\rm
obs}}$ (with exceptions for sources very close to merging) and thus the SNR of
a typical source in a 10-year LISA mission is approximately $1.58$
(=$\sqrt{10/4}$) times stronger than in a 4-year mission.

We use LEGWORK (Wagg et al., 2021) to calculate the signal-to-noise ratio for
each binary and the package ensures that enough harmonics are computed for
each binary such that the error on the gravitational-wave luminosity remains
below 1%.

####  2.3.3 Detection rate calculation

For each physics variation model and DCO type, we first convert the COMPAS
simulation results into a total number of DCOs in the Milky Way, $N_{\rm
DCO}$. We do this by integrating the full mass and period distributions and
stars and normalising to the total Milky Way mass. For more details see
Appendix B.

We then determine the fraction of binaries that are detectable in each Milky
Way instance by summing the adaptive importance sampling weights of the
binaries that have an SNR greater than 7, and dividing by the total weights in
the simulation. We multiply this fraction by $N_{\rm DCO}$ to find a detection
rate (which we write as a total number of detections per LISA mission)

$N_{\rm detect}=\frac{\sum_{i=0}^{N_{\rm MW}}w_{i}\phi(i)}{\sum_{i=0}^{N_{\rm
MW}}w_{i}}N_{\rm DCO},$ (13)

where $\phi(i)=1$ if a binary is detectable and $0$ otherwise. We calculate
the detection rate by Monte Carlo sampling 2500 Milky Way instances (each
containing 200,000 DCOs) for each DCO type and every physics variation in
order to obtain values for the uncertainty on the expected detection rate.

##  3 Results I - Predictions for LISA sources

In this section we present our predictions for the population of detectable
LISA sources for our fiducial model. In total we expect, on average, $124$
detections in a 4-year LISA mission, of which $74{}$, $42{}$ and $8{}$ are
BHBHs, BHNSs and NSNSs respectively, based on our fiducial simulations. In the
remainder of this section, we discuss where the sources are expected relative
to LISA’s sensitivity curve (Sec. 3.1), their properties (Sec. 3.2), their
locations in the Milky Way (Sec. 3.3), their formation channels (Sec. 3.4) and
finally we discuss the expected SNR and how accurately we expect that the
parameters can be measured (Sec. 3.5). Note that all results shown in this
section are based on our fiducial simulations. A discussion of the impact of
variations in the physics assumptions is provided in Sec. 4.

###  3.1 The LISA sensitivity curve and the population of detectable sources

Figure 2: The LISA sensitivity curve is shown together with the density
distribution of the characteristic strain and the dominant frequency for all
detectable sources in our simulations (top) and separated by type (bottom).
Contours show the percentage of the population enclosed. The remaining 2% of
the population is shown as dots with a size that scales with the sampling
weight. To guide the interpretation, we show reference lines marking where a
circular binary would reside for a given distance (diagonal line) and
remaining inspiral time (vertical lines), assuming an average chirp mass
$\left\langle\mathcal{M}_{c}\right\rangle$, orientation and sky location. The
coloured lines in the bottom panels show a contour that encloses 90% of the
population that is circular. LISA verification binaries are overplotted in the
top panel (star symbols). See also Fig. 17 and Sec. 3.1 for a discussion.  .

We show the expected LISA sensitivity curve (Robson et al., 2019) in Fig. 2,
which includes the confusion noise arising from the Galactic WDWD population,
and overplot our predictions for the distribution of detectable sources.
Eccentric systems emit gravitational waves in multiple harmonic frequencies
($nf_{\rm orb}$, with $n=2,3,\dots$). We choose to plot them at the
$x$-coordinate that corresponds to the frequency of the harmonic that
individually accumulates the largest SNR. For circular systems, the
$x$-coordinate simply corresponds to $2f_{\rm orb}$. The $y$-coordinate
indicates the strength of the signal (or to be more precise, the amplitude
spectral density, ASD), including the contribution from all harmonics.

For reference, we show dotted lines to indicate where a hypothetical binary
system would reside assuming a given distance from earth (diagonal lines) and
a fixed remaining inspiral time (vertical lines). For each line we assume the
binary is circular and has a chirp mass equal to the average of the sample
($\left\langle\mathcal{M}_{c}\right\rangle$, annotated in each panel). We also
overplot the LISA verification binaries (star symbols, Kupfer et al. 2018).

We observe several features in Fig. 2 that are worth discussing and
explaining. We note that some of these have also been described in earlier
studies (see Sec. 6). Firstly, we note the empty band that separates the LISA
sensitivity curve and the detectable population. This is the result of the
criteria for detection where we require $\rm{SNR}>7$.

We further note that the detectable population is concentrated on the left
side of LISA’s sensitivity window. The peak is located at a frequency of about
$0.4\,\mathrm{mHz}$, which is much lower than the frequency where LISA will be
most sensitive (about $10\,\mathrm{mHz}$). This can be understood from the
acceleration of inspiraling DCOs as they evolve towards higher frequencies.
DCOs typically form with wide orbits (low frequencies) that would not be
detectable yet. Their orbits shrink as they lose angular momentum in the form
of gravitational waves leading to an increase of their orbital frequency until
they become detectable. These systems are increasingly rare because they
evolve faster and faster towards higher frequency as the inspiral accelerates,
even though the signal emitted by a more compact binary is stronger (Peters,
1964). The vertical grid lines show these rapidly decreasing inspiral times at
increasing frequencies. Most of the population is thus expected to reside at
low frequencies.

In the lower three panels, we show the contributions of the different types of
sources. Comparing them, one can observe the shift in the frequency at which
the peak is located, at $0.3\,\mathrm{mHz}$, $0.7\,\mathrm{mHz}$ and
$1\,\mathrm{mHz}$ for BHBH, BHNS and NSNS systems respectively. This is a
result of the difference in chirp mass. A higher mass system can emit at lower
frequency and still produce the minimum SNR needed for detection. We note that
this effect can be used to distinguish the heavier DCOs that we discuss in
this work from WDWD systems, at least probabilistically (see Sec. 5.1.1). In
the same way, this also explains the offset in frequency between the
population of sources we predict and the LISA verification binaries.

Inspecting the dotted reference lines, we note that the peak of the density
distribution of observable sources coincides with the location expected for
circular systems at $8\,\mathrm{kpc}$, which is the distance to the centre of
the Milky Way. As can be seen best in the lower panels, the reference lines
for $0.1$ and $30\,\mathrm{kpc}$ enclose the majority of systems, as expected
given the dimensions of the Milky Way.

There is a distinct subpopulation of binaries, most clearly visible in the
lower panels as an offshoot that extends downwards to ASDs of
$10^{-19}\,\mathrm{Hz}^{-1/2}$, especially around $2\,\mathrm{mHz}$. This
offshoot is almost uniquely composed of eccentric binaries, as can be seen in
Fig. 17, which shows a similar figure but colouring individual systems by
their eccentricity. This can also be seen, albeit indirectly, from the contour
lines shown in the bottom panels of Fig. 2, which encompass 90% of the
circular sources in each population. This contour does not include the
offshoot. We conclude that eccentric sources occupy a very different region in
this diagram.

###  3.2 Properties of the detectable systems

In Fig. 3, we show the 1D distribution of several individual parameters of the
population of detectable binaries together with the 1- and 2-$\sigma$
uncertainties obtained via bootstrapping. These uncertainties represent the
fluctuations in our results over different random instances of the Milky Way.
The distributions shown here are approximated by kernel density estimators,
corrected for edge effects by mirroring the sample (Schuster, 1985).

##### Orbital Frequency

The orbital frequency distributions for BHBHs, BHNSs and NSNSs (Fig. 3a) peak
at progressively increasing frequencies. As mentioned in Sec. 3.1, this is
because a higher mass DCO at the same distance and eccentricity requires a
lower frequency to produce the same signal-to-noise ratio and thus be
detected. The distributions appear nearly symmetric, but closer inspections
shows that the left hand side is more populated, which can be seen most
clearly in the curve for the BHBHs. This is due to the contribution of highly
eccentric binaries, which are most abundant in the BHBH population. These
systems are still detectable by LISA, despite their low orbital frequency, as
the high eccentricity means that the majority of the GW signal is emitted at
higher harmonics, where LISA is more sensitive.

Figure 3: Properties of detectable systems for a 4-year LISA mission in our
fiducial model. Each panel shows a kernel density estimator for a single
property, coloured by DCO type. Each curve has been individually normalised.
The shaded areas show the 1- and 2-$\sigma$ sampling uncertainties (obtained
via bootstrapping). The dotted lines in panel b show the individual primary
and secondary mass distributions. The dotted line in panel e shows the star
formation history we have assumed in our Milky Way model. See Sec. 3.2 for a
discussion.  .

##### Black Hole Mass

In Fig. 3b, we show the distribution masses of individual black holes that are
part of BHNS and BHBH systems. The distribution for BHBH shows a bimodality.
This results from the two contributions of the more and less massive BHs in
BHBH systems, which peak at around $8\,\mathrm{M_{\odot}}$ and
$3.5\,\mathrm{M_{\odot}}$ respectively as shown by the dotted lines (see also
our discussion of the mass ratios below).

For both the BHBHs and BHNSs, we see that the black hole mass distribution
favours low masses. About $90\%$ of BHs have masses below
$11\,\mathrm{M_{\odot}}$, in our fiducial simulations shown here. This is in
stark contrast with observations from ground-based GW detectors, where heavy
BHs with masses of about $30\,\mathrm{M_{\odot}}$ and higher have been common.
There are two main reasons for this discrepancy. First, the population of BHs
in the Milky Way (and, in particular, those detectable by LISA due to their
recent formation times, see Fig. 3e) primarily come from progenitors that
formed from high metallicity gas according to our simulations. Stellar winds
are stronger at high metallicity leading to increased mass loss. This affects
the mass of the most massive black holes that can be formed (Belczynski et
al., 2010). Secondly, ground-based detectors are strongly biased towards high
mass systems, since they can be detected out to larger distances and thus a
greater volume is probed. In contrast, LISA has no such bias and is, in
principle, sensitive to inspiraling BHBHs throughout the entire Milky Way
regardless of their mass, as long as we catch them when they are emitting in
the LISA band. For this reason, LISA is more likely to detect the more common
lower mass BHBHs.

We also note that our BH mass distribution extends down below
$5\,\mathrm{M_{\odot}}$ to $2.5\,\mathrm{M_{\odot}}$ which is our fiducial
maximum neutron star mass. The BHs in this simulation fill the so-called
“lower mass gap” marked as a grey band (Özel et al., 2010; Farr et al., 2011),
see also Shao & Li (2021) who recently also pointed this out. This prediction
is sensitive to adopted model for fallback during the SN explosion as we
discuss this further in Section 4.2.1.

##### Mass Ratio

The mass ratio distributions for each DCO type are very distinct from one
another, as can be seen in Fig. 3c. The majority of NSNSs have a mass ratio
close to unity, with $90\%{}$ of systems having $q>0.8$, where $q\equiv
m_{2}/m_{1}$. The reason for the concentration around equal masses is that
most NSs are formed either through electron-capture supernovae (ECSN) or from
low mass stars in our simulations. We assume a remnant mass for any NS formed
through ECSN of $1.26\,\mathrm{M_{\odot}}$ (see Sec. A.2). The remnant mass
prescription that we use assumes a fixed fallback mass for any star with a CO
core mass less than $2.5\,\mathrm{M_{\odot}}$, such that many NSs end up with
an identical mass of $1.278\,\mathrm{M_{\odot}}$ (see Fryer et al., 2012, Eq.
19). This means that many NSs are formed with equal masses and hence we see a
mass ratio distribution peaked around unity.

In contrast, only $8\%{}$ of detectable BHBHs are formed with $q>0.8$ and the
distribution peaks around $q=0.4$. We find that the strong stellar winds in
our (typically high-metallicity) progenitors are the reason behind this.

BHBHs with unequal masses typically come from progenitors that also had more
extreme mass ratios at birth ($90$ and $30\,\mathrm{M_{\odot}}$ are typical
for the progenitors of detectable BHBH systems in our simulations. The primary
in such systems experiences strong mass loss by winds before filling its Roche
lobe. This mostly happens during its early hydrogen-shell burning phase. The
wind mainly reduces the mass of the envelope, but does not have a very strong
effect on the core. By the time the primary fills its Roche lobe, it has
become less massive and the mass ratio is closer to one. This favours stable
mass transfer. The massive core of the primary star typically becomes the more
massive BH. Accretion on the secondary star is limited and the secondary
eventually provides the less massive black hole.

At the same time, stellar wind mass loss disfavours the formation of black
holes with comparable masses. Such systems would have originated from
progenitors that also started with comparable masses. The rather massive
secondaries in these systems (especially after they accreted from the primary)
experience very strong stellar wind mass loss due to LBV-like eruptions. This
limits the amount by which they can expand (e.g. van Son et al., 2021). At the
same time, wind mass loss (and also stable and more conservative mas transfer
in these systems) lead to widening of the orbit. Both effects, the reduced
expansion and increased widening of the orbit, tend to prevent the secondaries
from being able to fill their Roche lobe. This thus limits the number of
systems that experience the reverse mass transfer or common-envelope phase
needed shrink the orbits and make them tight enough to be detected as
gravitational-wave sources.

We find that detectable BHNSs have even more unequal mass ratios. Moreover,
the mass ratio distribution is bimodal, where the two peaks arise from two
distinct formation scenarios. Around two thirds of detectable BHNSs experience
at least one common-envelope event, whilst the last third are formed through
only stable mass transfer. The first peak at $q=0.18$ is from systems that
experience at least one common-envelope phase and occurs at the expected mass
ratio, which approximately follows the mean BH mass
(${\sim}6.5\,\mathrm{M_{\odot}}$) and NS mass
(${\sim}1.2\,\mathrm{M_{\odot}}$). Yet we also see a second peak at higher
mass ratios around $q=0.34$, which arises from the fraction of the population
that underwent only stable mass transfer phases. The stability of mass
transfer depends on the mass ratio, preferentially forming systems with more
equal masses, i.e. at higher $q$.

##### Eccentricity

In Fig. 3d we show the eccentricity distributions. We find that most systems
(73%) will have eccentricities larger than $0.01$ during the LISA mission,
which should in principle be detectable according to Nishizawa et al. (2016).
This means that we will potentially be able to use eccentricity to distinguish
these sources from WDWDs, which are expected to have little to no eccentricity
(see Sec. 5.1.1). We note that several previous studies assumed all systems to
be circular when calculating the detection rates (e.g. Liu & Zhang, 2014;
Lamberts et al., 2018; Sesana et al., 2020). We discuss the impact of this
assumption in Section 6.

For systems with eccentricities higher than $e\gtrsim 0.3$, most gravitational
wave energy is emitted in higher harmonics. Such systems are more rare, but we
find them to be significant among the BHBH population, where they account for
21% of systems. Detectable BHBHs in our simulation (and, in particular, those
that are eccentric) are primarily systems formed through the stable mass
transfer channel (see Fig. 15). These systems are still relatively wide
(compared to those formed through the CE channel) immediately prior to
formation of the second BH, which makes them more easily affected by kicks. If
the kick is oriented roughly in opposite direction to the orbital motion and
has a velocity that is of similar magnitude as the orbital velocity, it will
lead to the formation of a highly eccentric system.

It is rare to get such a “lucky kick”, but there are a few effects that favour
this for BHBHs. The kicks of BHs are reduced by fallback and they are thus
less likely to disrupt the system. Moreover, because of their higher masses,
BHBHs can be observed already at lower orbital frequencies. This means that
they have not had as much time to circularise and so still have significant
eccentricity by the time of the LISA mission. Finally, LISA favours the
detection of eccentric systems, if all other properties are held fixed. This
is because the gravitational-wave emission is stronger (Eq. C2) and the energy
is emitted at higher frequencies (Peters & Mathews, 1963, Eq. 20) where LISA
is more sensitive.

The lower abundance of highly eccentric systems among the NSNS and BHNS
systems may seem counter-intuitive since neutron stars are lower mass and
would be more strongly affected by natal kicks, which one may expect to lead
to more eccentric systems. However, the majority of NSs in our simulations are
formed through ECSN and USSN and for these types of supernovae we draw from a
Maxwellian with $\sigma_{\rm rms}^{\rm 1D}=30\,\mathrm{km}\,\mathrm{s^{-1}}$.
Thus the kicks received by NSs in our simulations are often much smaller than
for BHs.

##### Time since formation

In Fig. 3e we show how long ago the LISA detectable DCOs formed. Star
formation was highest at early times $6-12\,\mathrm{Gyr}$ ago, after which it
declined. In contrast, the LISA detectable DCOs primarily formed in the
relatively recent history of our Milky Way, about $2\,\mathrm{Gyr}$ ago. This
reflects the fact that binaries in our simulation typically take about a Gyr
to merge.

When comparing the distribution of formation times for the three different DCO
types we see that NSNSs are most strongly concentrated at recent times,
followed by BHBHs and then BHNSs. To understand this it is helpful to consider
how the inspiral time scales with various parameters (Peters, 1964)

$t_{\rm
inspiral}\propto\frac{a^{4}}{(m_{1}+m_{2})^{3}}\cdot\frac{q}{(1+q)^{2}}\cdot(1-e^{2})^{-7/2}.$
(14)

The inspiral time depends most strongly on the separation at DCO formation,
$a$, and this is where the three types also differ most strongly (see Fig.
17). The detectable NSNS systems have the tightest orbits at DCO formation.
The median of the distribution of separations at DCO formation, $\langle
a_{\rm DCO}\rangle_{\rm med}$, relate as 8:3:1 for detectable BHBH:BHNS:NSNS
in our simulations. This results in increase of the inspiral time by a factor
of about 4000:80:1. The total masses affect the inspiral time to the third
power and this where the heavier BHBH systems are favoured. The median total
masses differ by ratios of 6:4:1 for detectable BHBH:BHNS:NSNS in our
simulations, impacting the inspiral times such that they are a factor of
200:60:1 shorter, partially counteracting the effect of the separations. The
term depending on the mass ratio $q$ only varies by about 30% for the mass
ratio ranges considered here and so is not of interest. The eccentricity term
is not of importance for mildly eccentric systems, $f(e_{\rm DCO}\leq 0.3)\leq
1.4$ but of large importance for the very eccentric $f(e_{\rm DCO}\geq
0.9)\geq 300$. The fraction of highly eccentric systems with $e_{\rm DCO}>0.9$
is 33%, 16%, and 8% of for BHBH, BHNS and NSNS respectively, see also Fig. 17.

We conclude that the shorter median separations at DCO formation are the main
reason why NSNS are most strongly peaked at short delay times. They are
followed by BHBHs rather than BHNSs due to the high masses and substantial
eccentricities of BHBHs.

##### Time until merger

Fig. 3f shows the remaining time until merger for each of the DCO types at the
start of LISA mission. The distributions are strikingly similar and peak with
merger times of around a Myr.

The merger time is a function of the mass, frequency and eccentricity of the
sources, such that more massive, higher frequency and more eccentric sources
merge faster (Peters, 1964, Eq. 5.14). So, despite the fact that each DCO type
often has higher values in any one of these properties, the convolution of all
three tends to negate the differences. For example, NSNSs have the highest
orbital frequencies and are mildly eccentric whilst BHNSs have moderate
orbital frequencies and are more circular. However, BHNSs are more massive in
general and so the overall merger times are distributed very similarly for
both DCO types.

###  3.3 Distribution in the Milky Way

For our Milky way model we consider three different components, a
low-$[\alpha/\mathrm{Fe}]$ (“thin”) disc, an older high-$[\alpha/\mathrm{Fe}]$
(“thick”) disc and a bulge/bar (see Sec. 2.2). In Table 1 we summarise the
number of detections originating from each of these components. Despite the
fact that only $42.5\%$ of systems are formed in the
low-$[\alpha/\mathrm{Fe}]$ disc, we find that $86\%$ of the detections
originate from this component. This is because most detectable systems were
formed relatively recently (see Fig. 3e) and so the
high-$[\alpha/\mathrm{Fe}]$ disc and bulge are effectively too old to
contribute many detectable systems. Nevertheless, we do find a significant
fraction of detectable systems originate in the high-$[\alpha/\mathrm{Fe}]$
disc and bulge, indicating that it is still important to include these
components, as was ignored in some earlier works (see Sec. 6).

Component | Formation | Detectable
---|---|---
All | All | BHBH | BHNS | NSNS
Low-$\alpha$ disc | 42.5% | 86% | 89% | 82% | 85%
High-$\alpha$ disc | 42.5% | 6% | 5% | 8% | 12%
Bulge/bar | 15.0% | 8% | 6% | 10% | 3%
Table 1: Percentage of systems in each Galactic component.

Figure 4:  Top: A face-on view of the Galactic density distribution for
detectable DCOs. We show the density distribution for the top 95% of the
sources, the rest are indicated by scatter points whose sizes correspond to
their sampling weights. The coloured lines show the 75% contour for each of
the individual DCO types. The large cross passes through (0, 0) and helps to
highlight the bias towards the position of the sun, which is indicated by the
$\odot$. Bottom: As Fig. 3, but for the luminosity distance.  .

In the top panel of Fig. 4, we show the density distribution for detectable
DCOs in the galaxy. We see that most detectable sources are concentrated
towards the Galactic centre, with a strong bias towards sources that are on
our side of the Milky Way in the vicinity of the solar system (indicated with
the $\odot$ symbol). In principle, systems are detectable out to large
distances of about $20\,\mathrm{kpc}$ and more, although they become
increasingly rare, as can be seen from the 95% contour.

The differences between different DCO types can be seen more clearly in the
bottom panel of Fig. 4 where we show the distribution of the distances, $D$,
from earth to the detectable systems. Each distribution peaks around
$8\,\mathrm{kpc}$, which is the distance to the centre of the Milky Way. The
distribution for BHBH and BHNS systems follow a very similar shape, favouring
the detection of systems with distance $<8\,\mathrm{kpc}$, but with a tail
extending out to about $20\,\mathrm{kpc}$.

The distribution for NSNS stands out by being flatter, making them more common
nearby and, surprisingly, also at larger distances compared to the stellar
density. This may seem counter intuitive as one might naively expect the less
massive NSNS systems would not be observable out to larger distances than the
more massive BHNS and NSNS systems. To understand the differences we need to
consider not only the mass distribution of binaries, but also their
eccentricity and frequency distributions. Since, each parameter contributes to
the calculation of the SNR (and thus affects the maximum distance at which
systems can be detected).

The reason that the NSNSs are favoured at higher distances is that the NSNS
population has the highest fraction of “mildly” eccentric systems
($0.01<e<0.03$). In contrast, the BHNS population has a much higher fraction
of effectively circular systems ($e<0.01$), which emit weaker gravitational
waves compared to equivalent eccentric systems. Therefore, despite their
typically higher masses, the distance at which a BHNS source is detectable is
generally lower than for the mildly eccentric NSNS.

Conversely, the BHBH population has the highest fraction of highly eccentric
systems ($e>0.3$). Although one may naively expect that this would result in
stronger signals (and so further distances), for a system to have these high
eccentricities in LISA, it must still be early in its evolution (otherwise it
would have circularised) and thus have a low orbital frequency. The result of
this is that highest eccentricity systems tend to have lower SNRs and so
cannot be detected at large distances.

Overall we see that the eccentricity distribution of NSNSs occupies a “sweet
spot” where the gravitational wave power is increased compared to circular
systems, but it isn’t too high that the frequency is significantly impacted.
This means that NSNSs can be seen out to the largest distances of the three
DCO types.

###  3.4 Formation channels

In our fiducial model, approximately two thirds of detectable BHBHs are formed
through the ‘only stable mass transfer’ channel, whilst the remaining third
are primarily formed through the ‘classic’ CE channel. Detectable BHNSs follow
an inverse pattern, such that around two thirds are formed through the classic
channel and the rest are mainly formed through only stable mass transfer.

In contrast, detectable NSNSs are very rarely formed through only stable mass
transfer. Approximately half of systems are formed through the ‘classic’
channel and the rest are formed through a double-core common-envelope event
(Brown, 1995) where both progenitors evolve on a similar timescale and
initiate a double-core common-envelope event whilst they are on the giant
branch. All detectable DCOs show a small fraction of systems are formed
through a channel that does not fit into the other categories and hence are
labelled ‘other’. These systems tend to be formed through ‘lucky’ supernova
kicks that happen to shrink the binary significantly by chance.

The fraction of detectable DCOs that are formed through different formation
channels in the other model variations are shown Fig. 15, where the first
column in the plot corresponds to the fiducial model that we described above.

###  3.5 How accurately will we be able to infer the parameters of detected
systems?

So far we have discussed the properties of the all detectable sources.
However, only for a subset of systems we expect to get high enough SNR to
obtain accurate and useful measurements of these parameters. Below we discuss
the expected SNR distribution and the typical uncertainties expected for the
most relevant parameters, namely, the chirp mass and sky localisation.

##### Signal-to-noise ratio

Figure 5: Cumulative number of LISA detections with a given signal-to-noise
ratio. Colours indicate LISA mission length and shading shows 1- and
2-$\sigma$ uncertainties (obtained via bootstrapping).  .

In Fig. 5 we show the cumulative number of detections with a given SNR.
Although many a large fraction of sources have SNRs around our assumed
detection threshold of 7, many systems are detected with very high SNRs. We
find that on average for a 4(10)-year LISA mission, of the 124(202)
detections, 85(138), 16(27) and 9(14) systems have SNRs greater than 10, 50
and 100, respectively. These high SNR systems are typically, but not only, the
more massive BHBH systems.

##### Chirp mass

The chirp mass is important for identifying the type of the source of a
detected GW signal. The uncertainty of the chirp mass depends on the
uncertainty in the measured orbital frequency, the time derivative of the
orbital frequency and the eccentricity as detailed in Appendix C.

We find that for a 4(10)-year LISA mission, approximately 41(105) detections
have measurable chirp masses ($\Delta\mathcal{M}_{c}/\mathcal{M}_{c}<1$,
indicated by the dark shaded region in Fig. 6) whilst 13(31) have chirp mass
uncertainties below 10% ($\Delta\mathcal{M}_{c}/\mathcal{M}_{c}<0.1$,
indicated by the light shaded region in Fig. 6). This uncertainty is generally
dominated by the uncertainty on the time derivative of the frequency, since
most of the binaries are too early in their inspiral for LISA to measure a
strong chirp. Note from Fig. 6 that increasing the mission length
significantly increases the number of detections for which the chirp mass
uncertainty is below 100%. The total number of detections only scales as
$\sqrt{T_{\rm obs}}$, yet we find that the the total number of detections with
a chirp mass uncertainty is below 100% and 10% both scale approximately as
$T_{\rm obs}$.

Figure 6: As Fig. 5, but for the chirp mass uncertainty. The shaded areas
indicate regions with more than 10% and 100% uncertainty.  .

##### Sky localisation

An accurate sky localisation will be essential to possibly identify
electromagnetic counterparts or distinguish sources that come from different
components of our Milky Way.

We quantify the sky localisation of a source by estimating the angular
resolution for the detectable sources. Since all potential sources are
effectively stationary on the timescale of the LISA mission, we can follow
Mandel et al. (2018) and use the timing accuracy of LISA and the effective
detector baseline to calculate the angular resolution, $\sigma_{\theta}$, as

$\sigma_{\theta}=16.6^{\circ}\left(\frac{7}{\rho}\right)\left(\frac{5\times
10^{-4}\,\mathrm{Hz}}{f_{\rm
dom}}\right)\left(\frac{2\,\mathrm{AU}}{L}\right),$ (15)

where $L$ is the effective detector baseline, which for LISA is 2 AU in the
frame of the solar system, since it will orbit the Sun.

We plot the distribution of expected angular resolutions in Fig. 7. We see
that, for a 4(10)-year LISA mission, approximately 82(123) sources can be
resolved to an angular resolution better than 10 degrees and only 16(23)
better than 1 degree.

For comparison, the size of a pencil beam for a $15\,\mathrm{m}$ diameter SKA
dish observing at 1.4 GHz is roughly 0.67 square degrees (Smits et al., 2009),
corresponding to an angular resolution of
$\sigma_{\theta}=\sqrt{(0.67/\pi)}=0.46^{\circ}$ (similar to the angular size
of the moon). We will further discuss the prospects of matching LISA
detections to radio pulsars with SKA in Sec. 5.2.

Figure 7: As Fig. 5, but for the angular resolution.  .

##  4 Results II - Impact of physics assumptions

In this section we explore the effect of varying the uncertain assumptions
governing the evolution of binary system and the formation of compact objects.
For this we use the population synthesis simulations and model variations
presented in Broekgaarden et al. (2021) and Broekgaarden et al. 2021b (in
prep.).

We first discuss the robustness of our predictions for the number of
detectable systems (Section 4.1). We then discuss examples of how the
observable properties of the detectable systems, such as the distribution of
masses and eccentricities, are affected and can potentially be used to probe
the physics of double compact object formation (Section 4.2).

###  4.1 Detection rates

We predict approximately 30-370 detections in a 4-year LISA mission, across
all our simulations for varying physics assumptions. This increases to about
50-550 for a 10-year LISA mission. Although the number of detections per type
can vary by about 2 orders of magnitude, we find that the total detection rate
is fairly robust, among the variations we have considered (see Table 3).

In Fig. 8, we show the expected number of LISA detections based on our
simulations considering variations in the physical assumptions. We show the
expected number of detections for BHBH, BHNS and NSNS systems in the top,
middle and bottom panel respectively. All the rates and their uncertainties
plotted in this figure are also provided in Table 3. In the sections that
follow, we briefly explain the variations considered and discuss the most
prominent trends.

####  4.1.1 Efficiency of mass transfer

The efficiency of mass transfer, i.e. the fraction of mass lost by the donor
through Roche-lobe overflow that is accreted by the companion, is poorly
constrained and is considered as one of the main uncertainties in binary
evolution (e.g. de Mink et al., 2007). In our fiducial model A, we use a
prescription in which the accretion rate onto stellar companions is regulated
by their thermal timescale, i.e. the timescale on which a star can react to
changes and restore thermal equilibrium (see e.g. Schneider et al., 2015).

In models B-D, we instead adopt a fixed value for the mass transfer
efficiency, $\beta$, from $\beta=0.25$ up to 0.75, in cases of stable mass
transfer onto a stellar companion. For accretion onto NS and BH we still
assume that their accretion is limited to the Eddington rate.

Nearly all systems that can be detected form through channels where the very
first interaction is stable mass transfer. Generally, higher mass transfer
efficiencies lead to higher masses for the accreting stars, but also leads to
wider orbits (Soberman et al., 1997; van Son et al., 2020). Changing $\beta$
thus already affects the masses and orbital separation after the first
interaction phase, which in turn changes the starting conditions and outcome
of all subsequent interactions phases. This makes it complicated to fully
understand the impact of varying $\beta$ in simple terms, but we can
distinguish two main patterns for higher and lower mass systems respectively.

For the most massive progenitors, increasing $\beta$ leads to secondary stars
that are so massive and luminous that they experience strong wind mass loss.
This leads to further widening of the orbit. In addition the more massive
secondaries may not be able to fully expand as mass loss may prematurely
remove their hydrogen envelope. Both of these effects tend to prevent the most
massive secondaries from filling their Roche lobe. This means that they cannot
initiate the reverse interaction needed to shrink the binary system and
eventually produce a detectable double compact object. We indeed see that
increasing $\beta$ leads to a decrease of the expected number of detections
for BHBHs and BHNSs, which originate from the most massive progenitors.

For lower mass progenitors, which primarily produce NSNS systems, we find the
opposite: increasing $\beta$ leads to an increase in the number of detectable
NSNS systems. This is in part because the changes in secondary mass and the
orbital widening, affect the number of systems for which the reverse
interaction successfully ejects the envelope and shrinks the orbit.
Furthermore, the increased mass of the secondary stars allows stars that would
have otherwise ended their life as a WD to instead become massive enough to
form a NS (e.g. Zapartas et al., 2017). The same effect also allows some NS
progenitors to become massive enough to become BH progenitors, which partially
cancels the extra progenitors that would have originally been destined to
become WDs.

Figure 8: The number of expected detections in the LISA mission for different
DCO types and model variations. Error bars show the 1- (solid) and 2-$\sigma$
(dotted) Poisson uncertainties. An arrow indicates that the error bar extends
to zero. The left axis and grid lines show the number of detections in a
4-year LISA mission and the right axis shows an approximation of the number of
detections in a 10-year mission (we scale the axis by $\sqrt{T_{\rm obs}}$,
see Table 3 for exact rates). Each model is described in further detail in
Table 2 and details of the fiducial assumptions are in Section A.2. See Sec.
4.1 for a discussion.  .

####  4.1.2 Common-envelope evolution

The common-envelope phase constitutes a highly uncertain phase in the
evolution of interacting binary systems (e.g. Ivanova et al., 2013). The
uncertainties concerns the conditions required for the onset of a common-
envelope phase and, if a common-envelope phase occurs, what the outcome is.
Rapid population synthesis simulations such as ours approximate both questions
in a crude way. We therefore consider several model variations.

##### The efficiency parameter $\alpha_{\rm CE}$

To estimate the outcome of a CE phase, we use a simple consideration of the
binding energy and orbital energy (Webbink, 1984; de Kool, 1990). Our fiducial
model assumes a common-envelope efficiency parameter $\alpha_{\rm CE}=1$ which
can be interpreted as the case where all the energy liberated by shrinking the
orbit is used in an optimal way to unbind the envelope. There have been many
attempts to constrain this parameter using observations and more recently also
using 3D simulations (e.g. De Marco et al., 2011; Law-Smith et al., 2020; Lau
et al., 2021), but no consistent picture has emerged. Therefore we consider
large variations in this parameter.

In models G-J we alter the common-envelope efficiency parameter $\alpha_{\rm
CE}$ to $0.1,0.5,2.0$ and $10.0$ respectively. Values smaller than 1 may
represent cases where not all energy is used efficiently to unbind the
envelope, for example when part of the energy escapes in the form of radiation
or if part is used to impart additional kinetic energy in the ejecta (e.g.
Ivanova et al., 2013; Nandez & Ivanova, 2016). Values larger than 1 may
represent cases where additional energy sources can be tapped into, such as
for example jets powered by accretion (e.g. Schreier et al., 2021). The
variations can also be seen as a way to cover uncertainties in estimates for
the binding energy itself.

Increasing $\alpha_{\rm CE}$ makes common-envelope ejection more efficient or,
in other words, less orbital shrinkage is needed to successfully eject the
envelope. This has two consequences. (1) A larger fraction of systems avoids
merging during the CE phase. This increases the overall number of DCOs. (2)
The systems that survive the CE phase are wider, possibly too wide to become
detectable as gravitational wave sources (e.g. Klencki et al., 2021). So,
while the first consequence favours the formation of DCOs, the second
consequence disfavours the formation of DCOs that are tight enough to be
detected.

These two opposing effect result in a fine tuning situation. Only a very small
subset of progenitor systems have the right orbital parameters prior to the CE
phase to successfully produce detectable systems. Changing $\alpha_{\rm CE}$
moves and changes this window in the parameter space that successfully leads
to the formation of detectable systems. How the number of the detectable
systems changes depends on whether the relevant parameter space grows or
shrinks and on how well the relevant part of the parameters space is
populated. Fully unravelling these effects and how they interplay with the
assumed star formation history is beyond our scope (and possibly not even of
large relevance given the severe simplifications). We will limit the further
discussions to simply stating the trends we observe.

We find that the BHBH rate peaks for $\alpha_{\rm CE}=1$ (model A) and reduces
whether we increase or decrease $\alpha_{\rm CE}$. The BHNS rate follows the
same pattern as BHBHs although the value of $\alpha_{\rm CE}$ which maximises
BHNSs seems to be between 1 and 2. In contrast, for NSNSs we find that
increasing $\alpha_{\rm CE}$ (models I and J) results in significantly higher
rates.

##### The “optimistic” CE treatment

We further explore a model variation introduced by Belczynski et al. (2007)
often referred to as the “optimistic” CE scenario. This variation (model K)
relaxes our restriction that donor stars that are on the Hertzsprung cannot
survive common-envelope events.

In agreement with other studies, we find that this treatment leads to a
significant increase in the formation rate of BHBHs, by a factor of two. This
is because the progenitors expand significantly during the Hertzsprung gap
phase in our simulations. In our fiducial simulation, all progenitors that
initiate unstable interaction during this phase would end as stellar mergers,
while in this variation they will survive. The progenitors of BHNSs and NSNSs
are less strongly affected with an increase about 30% increase.

##### Case BB mass transfer

In models E′ and F we consider uncertainties in case BB mass transfer. This is
a phase of mass transfer where the donor star has already lost its hydrogen
envelope in a prior interaction, but fills its Roche lobe again as it expands
during helium shell burning phase (e.g. Dewi et al., 2002; Tauris et al.,
2015, 2017). This is of particular interest for the formation of NSs, as their
lower mass progenitors are swelling the most during this phase (e.g. Laplace
et al., 2020, and references therein). Population synthesis studies find that
nearly all NSNS systems form through a phase of case BB mass transfer (Vigna-
Gómez et al., 2018).

In model E′ we enforce that case BB mass transfer is always unstable, such
that it always leads to a CE. Note that this is slightly different from model
E described in Broekgaarden et al. (2021). In their work the pessimistic
approach to CE evolution is implemented such that all HG stars are excluded
including helium stars in the helium shell burning phase. This, in combination
with the assumption that case BB mass transfer is always unstable, effectively
leads to the exclusion of all systems that originate through this channel. We
are interested in NSNS systems, which are frequently formed through this
channel in our simulations. Therefore, we adapted this model to only exclude
H-rich HG donors, but allow systems with donors that are helium stars in the
helium shell burning phase to survive a CE phase.

As expected, in model E′, we find that case BB systems form through a common-
envelope phase rather than only stable mass transfer (see Fig. 15). This model
gives the lowest NSNS rate of all of our variations, which is a factor of 3
lower than our fiducial rate. We also find a reduction of the BHNS systems by
a factor 2. The BHBH systems are not significantly affected, as expected,
since case BB mass transfer does not play a role for high mass progenitors.

Finally, in model F, we again enforce that case BB mass transfer is always
unstable, but in combination with the optimistic treatment for CE (essentially
combining models E′ and K). This allows the systems that have HG donors for
common-envelopes (as well as those formed through case BB mass transfer) to
survive the CE phase. We find that this model leads to the highest predictions
for the detections among all variations that we have considered.

####  4.1.3 Supernovae and compact remnants

The formation of compact remnants and their associated natal kicks also
constitute important uncertainties.

In model L we consider the so-called “rapid” remnant mass function by Fryer et
al. (2015) as an alternative to the “delayed” description used in our fiducial
simulations. This affects the mass distribution (see Sec. 4.2.1), but the
effects on the number of detectable systems is modest for BHBH and BHNS, and
negligible for NSNS.

The same is true for the impact of changing the assumed maximum neutron star
mass, $m_{\rm NS,max}$ (M and N). Lowering $m_{\rm NS,max}$ increases the
number of detections involving BHs (since more stars form BHs instead of NSs)
and vice versa, but has no significant effect on the number of NSNS detections
since the vast majority are formed from low mass NSs.

We find that not implementing pair-instability supernovae (PISN) or
pulsational pair-instability supernovae (PPISN) in model O has no effect on
the number of detections with LISA. This is because the average metallicity of
the Milky Way is high enough such that no progenitor retains enough mass to
initiate a PISN or PPISN.

Decreasing the natal kicks for all core-collapse supernovae (models P-Q)
increases the detection rates for each DCO type, since lower kicks result in
fewer disrupted binaries and hence a more numerous detectable population. The
BHNS and NSNS systems are strongly affected whilst the impact on BHBH systems
is insignificant. The reason that BHBHs are relatively unaffected is that, in
our models, the natal kicks for BHs are scaled down with the amount of mass
that falls back. In the case of BHBHs, the black holes have very massive cores
and thus low kicks.

In model R we assume BHs form without any kick, while using our fiducial
assumption for the natal kicks of neutron stars. This increases the
predictions for BHNS by a factor 3 but the impact on BHBH systems is much
smaller for the same reason as for models P-Q. As expected, the NSNS
population is not affected.

####  4.1.4 Stellar winds

Mass loss in the form of stellar winds or eruptions is also a main
uncertainty. It affects by how much stars can grow in size and it affects the
final core masses. We consider variations in the mass loss by naked helium
(Wolf-Rayet like) stars and choose to vary the efficiency of these winds
between $0.1$ and $5$, to account for uncertainties in the derived rates (e.g.
Vink, 2017; Shenar et al., 2019; Hamann et al., 2019; Sander & Vink, 2020).

We find that a reduction of the wind mass loss has very little effect on our
predictions. This is the consequence of several effects that cancel each
other. Firstly, the decreased Wolf-Rayet like winds mean that the DCOs
(particularly those containing BHs) are generally more massive and so more
detectable in LISA. Secondly, one may expect that LISA sources in model S
would be higher frequency than in our fiducial model as decreased winds
generally result in tighter binaries. However, though this is the case at DCO
formation, we find that by the time the sources have evolved until they are
observable by the LISA mission, they have lower orbital frequencies than in
our fiducial model. This is because the reduced winds allow DCOs to be formed
at higher metallicity and, therefore, at more recent times. This means that
most DCOs do not evolve for as long before the LISA mission and so remain at
lower frequencies (wider separations) thus making them less detectable.

In addition, we find that NSNSs are more eccentric and BHBHs are less
eccentric than our fiducial model (with BHNSs relatively unchanged). The
increase in eccentricity for NSNSs comes from the same reason as the lower
frequency, more recent birth times mean that binaries have less time to
circularise. The same is not true for BHBHs as the more massive systems are
less affected by supernova kicks and so fewer high eccentricity systems are
formed. Overall, despite the large differences in the system properties, these
three effects in combination leave the detection rates relatively unchanged.

In model T we instead increase the efficiency of Wolf-Rayet winds by a factor
of 5. In this model the detection rate of BHBHs decreases by over a factor of
10 and BHNSs by over a factor of 2 whilst the NSNS rate is relatively
unchanged. Increasing the efficiency of WR winds widens the orbit and
decreases the final masses of DCOs. This means that some progenitors that
would have formed LISA sources under our fiducial assumptions would not have
enough mass to produce a DCO, or produce a NS instead of a BH. The effect is
strongest for the progenitors of black holes, since these are most strongly
affected by Wolf-Rayet winds. This effect is less pronounced in NSs since the
rate of mass loss for the progenitors of these systems is low enough that
changing by a factor of 5 still does not impact the final fate. Moreover, DCOs
that are formed tend to be less massive and therefore less detectable.

###  4.2 Properties of detectable systems

In this section, we consider how varying underlying physics assumptions
changes the properties of detectable systems. We focus on several key
differences across physics variations rather than showing the differences in
every model and thus this section is by no means exhaustive.

####  4.2.1 Using LISA to investigate the lower mass gap

Figure 9: Comparison of the component mass distribution of LISA detectable
DCOs when using the Fryer delayed (model A) and rapid (model L) remnant mass
prescriptions. Distributions are plotted in the same way as Fig. 3, except all
DCOs types are shown in one curve and each type is weighted by its detection
rate in the respective model.  .

In Fig. 9, we show the component mass distribution for all LISA detectable
DCOs (BHBHs, BHNSs and NSNSs) for two different remnant mass prescriptions.
The grey distribution uses the Fryer delayed remnant mass prescription (Fryer
et al., 2012), which is our fiducial assumption (model A). This prescription
produces compact objects in the lower mass gap ($2.5\,\mathrm{M_{\odot}}\leq
m\leq 5\,\mathrm{M_{\odot}}$) and indeed we find that, of the LISA detectable
DCOs, approximately 69% of BHBHs, 39% of BHNSs and 0% of NSNSs have at least
one component in the lower mass gap. Overall, weighting by the relative
detection rates, this gives that, in our fiducial model, 55% of our predicted
LISA DCO detections would have at least one component in the lower mass gap
when using this remnant mass prescription. This equates to approximately 69
systems being detected in the lower mass gap. Alternatively, the blue curve in
Fig. 9 shows the same distribution but for the rapid remnant mass prescription
(Fryer et al., 2012), which we use in model L. In this case, no compact
objects are formed (and therefore, detected) in the lower mass gap.

From the stark difference between these models, it is clear that it is
difficult at this point to say with any certainty what fraction of systems
LISA will detect in the lower mass gap given the highly uncertain formation
rate of systems in this mass range. Indeed, we find that the percentage of
detectable systems with at least one component in the lower mass gap varies
between approximately 30-70% (or, in terms of detections, from 15 to 156) for
the different model variations. However, it is important to highlight that if
DCOs are formed with components in the lower mass gap, LISA will be able to
detect them. And thus, LISA could be a useful instrument for providing
constraints on the existence or non-existence of a lower mass gap based on the
mass distribution of detected DCOs.

####  4.2.2 Effect of natal kicks on eccentricity distribution

Figure 10: As Fig. 3d, but for model Q. For comparison, we show the mean
distribution for the fiducial model (model A) as dashed lines.  .

In Fig. 10, we investigate how decreasing the magnitude of natal kicks from
core-collapse supernovae affects the eccentricity distribution of LISA
detectable DCOs. For reference, we show the mean fiducial distributions (model
A) as dashed lines (see Fig. 3d for full comparison). In the main curves, we
reduce the velocity dispersion for core-collapse supernovae from
$265\,\mathrm{km}\,\mathrm{s^{-1}}$ to $30\,\mathrm{km}\,\mathrm{s^{-1}}$
(model Q).

We find that the LISA detectable BHBHs are significantly less eccentric with
weaker kicks, such that the population above $e=0.2$ is nearly completely
eliminated. This is because BHBHs are often massive enough to withstand strong
natal kicks without disrupting and these kicks tend to impart significant
eccentricity. In model Q, very few systems are ever given such strong kicks
and thus very few BHBHs are detected with significant eccentricity.

Since BHNSs are less massive than BHBHs and have more unequal mass ratios,
they are more vulnerable to disruption during supernova kicks. BHNSs can only
withstand strong kicks when they are aimed in the correct direction and so
only a small ‘lucky’ fraction of the fiducial population is highly eccentric.
Therefore in model Q, although we see that the population of highly eccentric
BHNS systems is eliminated (similar to BHBHs), the peak of the distribution
actually shifts to higher eccentricity. This is because a larger fraction of
systems are given weaker kicks that BHNSs can withstand and these impart much
more moderate eccentricities.

Finally, we find that the NSNS distribution is relatively unchanged between
model A and Q. This is not surprising however since the majority of NSNSs are
formed through electron-capture supernovae and ultra-stripped supernovae and
for these types of supernovae we use $\sigma^{\rm 1D}_{\rm
rms}=30\,\mathrm{km}\,\mathrm{s^{-1}}$ already (see App. A.2) and thus there
is little difference between the models.

Overall, compared to our fiducial model, we find that decreasing supernova
natal kicks, though it strongly increases the number of detections (see Fig.
8), strongly decreases the fraction of highly eccentric systems that are
detected.

####  4.2.3 Effect of Wolf-Rayet winds on mass distribution

Figure 11: As Fig. 9, but instead varying Wolf-Rayet wind efficiency (models S
and T). The curve for $f_{\rm WR}=5$ has much higher uncertainties as there
are many fewer DCO systems formed in this model (the inverse reasoning also
explains the lower uncertainties for $f_{\rm WR}=0.1$).  .

In Fig. 11 we show the effect that changing the efficiency of Wolf-Rayet winds
has on the individual component mass distribution. Decreasing the Wolf-Rayet
wind efficiency allows the formation of more massive DCOs in the Milky Way
and, indeed, we see that the distribution extends to $25\,\mathrm{M_{\odot}}$
and relatively fewer detectable systems are formed at low masses. By contrast,
increasing the Wolf-Rayet efficiency by a factor of 5 strongly disfavours the
formation of systems at high masses and approximately 85% of detectable
systems have masses below $5\,\mathrm{M_{\odot}}$. These three distributions
are very distinct and so it is possible that the mass distribution of LISA
could help to constrain the efficiency of Wolf-Rayet winds.

##  5 Discussion

In this section we discuss the prospects of (and methods for) identifying LISA
sources (Sec. 5.1), the possibility of matching LISA signals to SKA detections
(Sec. 5.2), the main caveats for this study (Sec. 5.3) and the possible
contribution from other formation channels (Sec. 5.4). All predictions quoted
in each subsection are derived for the fiducial model (model A).

###  5.1 Identification of GW sources

It is important to note that, though we present predictions for the detection
rates of specific DCO types, the nature of the source may not be immediately
apparent from the gravitational wave signal. LISA can detect a variety of
sources, from exoplanets (e.g. Tamanini & Danielski, 2019) to common-envelopes
(e.g. Ginat et al., 2020; Renzo et al., 2021) that may cause confusion.
However, by far the most prominent will be the population of Galactic WDWDs
detectable with LISA, which will be several orders of magnitude larger than
the population of the more massive DCOs that we focus on in this paper (e.g.
Korol et al., 2017). It is therefore imperative that we consider how to
distinguish NS and BH binaries from this much more numerous population of
sources. In addition to distinguishing them from WDWDs, we must consider how
to discriminate between BHBHs, BHNSs and NSNSs themselves.

####  5.1.1 Distinguishing from WDWD population

The simplest way to check whether a source is a WDWD is to evaluate its chirp
mass. The mass of a non-rotating white dwarf cannot be larger than the
Chandrasekhar limit of $1.4\,\mathrm{M_{\odot}}$ (Chandrasekhar, 1931; Hamada
& Salpeter, 1961), so we can take the maximum chirp mass of a WDWD to be
${\sim}1.2\,\mathrm{M_{\odot}}$. Therefore, any DCO with a chirp mass that
satisfies $\mathcal{M}_{c}>1.2\,\mathrm{M_{\odot}}+\Delta\mathcal{M}_{c}$ must
not be a WDWD (where $\Delta\mathcal{M}_{c}$ is the error on the chirp mass,
estimated using Eq. C3). We find that for the detectable population of a
4(10)-year LISA mission, 24(38)% of BHBHs, 28(41)% of BHNSs and 4(5)% of NSNSs
satisfy this condition. This method is not particularly effective for NSNSs
since their average chirp mass, $1.17\,\mathrm{M_{\odot}}$, is below the
Chandrasekhar limit.

Another discriminator between WDWDs and other DCOs is eccentricity. WDWDs
formed in the disc are thought to be formed mainly through isolated binary
formation and have little to no eccentricity (e.g. Nelemans et al. 2001, see
however Dosopoulou & Kalogera 2016a, b; Gosnell et al. 2019). This is because
WDWDs formed through isolated binary evolution all experience a phase of mass
transfer or a common envelope, which typically circularises the binary (e.g.
Marsh et al., 2004). However, in contrast to the more massive DCOs that we
study, WDWDs do not experience strong natal kicks which we find to be the main
source of eccentricity. Therefore, if any system is detected with anything
other than one detectable harmonic, this suggests that the system is unlikely
to be a WDWD. We find that for a 4(10)-year LISA mission, 55(61)% of BHBHs,
27(29)% of BHNSs and 66(68)% of NSNSs are detected with multiple harmonics
(see also Sec. 3.2). Both the absolute percentage and the relative improvement
with an extended LISA mission is lower for the BHNSs with respect to other
DCOs as we find that these BHNSs are less eccentric on average (see Fig. 3d
and discussion in Sec. 3.2).

However, we should also consider that eccentric WDWDs could be formed through
dynamical formation in Milky Way globular clusters (e.g. Willems et al., 2007;
Kremer et al., 2018), or with third companions (e.g. Antonini et al., 2017).
This means that we cannot assume that eccentric binaries are not WDWDs unless
they are detected in the Galactic plane (though even then there is a chance
they were formed dynamically). We can use the sky localisation, scale height
of the disc and distance to the source to estimate what fraction of eccentric
sources can be localised to the Galactic plane. This condition can be written
as $\sigma_{\theta}<\arcsin(z_{\rm plane}/D_{L})$ or $D_{L}<z_{\rm plane}$,
where we set the height of the Galactic plane, $z_{\rm
plane}=0.95\,\mathrm{kpc}$, to the scale height of the high-$\alpha$ disc. We
apply this condition to find that the fraction of sources that are eccentric
and localised within the disc for a 4(10)-year LISA mission are 40(40)% for
BHBHs, 23(23)% for BHNSs and 59(59)% for NSNSs. Note that although the
fractions are the same for the 10-year mission, the absolute number of
detections is still greater.

Overall, combining these methods (chirp mass, eccentricity and sky
localisation) we find that for a 4(10)-year mission, LISA will detect at least
37(70) BHBHs, 18(38) BHNSs and 5(8) NSNSs that are distinguishable from the
WDWD population. Thus we will be able to confidently distinguish approximately
half of all detected sources from WDWDs. This increases to roughly 60% for a
10-year mission. We highlight that, though the overall number of LISA
detections in an extended mission only increases by a factor of $\sqrt{T_{\rm
obs}}$, the number of distinguishable detections increases by a greater factor
since each of the more numerous sources are better measured. This further
underlines the benefits of extending the LISA mission to 10 years.

####  5.1.2 Discriminating between BHBHs, BHNSs and NSNSs

The problem of discriminating between the BHBH, BHNS and NSNS populations can
be more difficult than distinguishing them from WDWDs. For NSNSs, we can
follow a similar method to the WDWDs (see Sec. 5.1.1) by applying our
knowledge of the maximum mass of a neutron star. Following our fiducial
assumption, we can take the maximum mass of a neutron star as
$2.5\,\mathrm{M_{\odot}}$ and thus the maximum chirp mass that a system can
attain without one of the components being a black hole is
$\mathcal{M}_{c}=2.2\,\mathrm{M_{\odot}}$. For a 4(10)-year LISA mission, the
fraction of systems that are above or below this limit (and thus must
respectively contain or not contain a BH component) by more than
$\Delta\mathcal{M}_{c}$ is 21(33)% for BHBHs, 18(24)% for BHNSs and 47(62)% of
NSNSs, which in terms of absolute detections is 16(39) for BHBHs, 7(17) for
BHNSs and 4(8) for NSNSs.

For separating the BHBH and BHNS population one could do so probabilistically
given the properties that are measured, particularly the orbital frequency,
mass ratio and eccentricity, since these distributions are fairly different
for the two DCO types (see Fig. 3). This method would pose a challenge,
however, as it would likely only indicate which type was more likely rather
than discriminate between them with strong evidence.

Another possible solution would be the existence of electromagnetic
counterparts to the gravitational wave signal. In Section 5.2 we consider the
possibility of detecting a pulsar within a BHNS or NSNS system. This could be
used to identify the type of the source.

###  5.2 Matching LISA detections to pulsars with the SKA

Since the vast majority of the LISA detectable population of DCOs will not
merge for many years, the main type of electromagnetic counterpart for this
population is pulsars. Therefore, for this section we focus only on BHNSs and
NSNSs since no BHBH system will contain a pulsar. The joint detection of a
binary pulsar with LISA and the Square Kilometre Array (SKA, Dewdney et al.
2009) would not only help to constrain the parameters of the binary, but also
enable investigation of other compact object physics. A pulsar(PSR)+BH can
provide stringent tests of theories of gravity, in particular the “No-hair
theorem” (Keane et al., 2015). Alternatively, an ultrarelativstic PSR+NS
system could be used to measure the neutron star equation of state up to an
order of magnitude more accurately than other proposed observational
constraints (Kyutoku et al., 2019; Thrane et al., 2020).

We estimate that on average, given the number of detectable pulsars and the
SKA sky area, each pulsar in the SKA occupies a region with an angular
resolution of $\sigma_{\theta}<1.3^{\circ}$ or $0.7^{\circ}$ for SKA-1 and
SKA-2 respectively (see Appendix E). Therefore, any DCOs containing NSs
localised by LISA with an angular resolution lower than these values can be
unambiguously matched to the radio signal in the SKA. By considering Fig. 7,
approximately $11$ and $6$ (for SKA-1 and SKA-2) DCOs will satisfy this
constraint.

If there is more than one pulsar in the region given by the LISA sky
localisation, one can compare the measured parameters of the system in LISA
and the SKA. Both the SKA and LISA will measure the orbital frequency to high
precision, as well as the time derivative of the frequency and chirp mass to a
lesser precision, of each of these systems. Therefore, one could perform a
targeted search with the SKA that checks the sky location given by LISA, only
looking for binary pulsars with orbital frequencies within the uncertainties.
If there was still more than one possible pulsar one could also check against
the chirp mass. In this way, we expect it will be possible to get a joint
detection between the SKA and LISA even when the sky area implied by the LISA
detection contains more than one pulsar.

In order to assess the efficacy of this method, we would need to know the
probability that two random binary pulsars would have orbital frequencies and
chirp masses close enough that one could not tell which pulsar matches the
LISA detection. This would require simulating the SKA population of pulsars
with a code such as PSRPOPPy (Bates et al., 2014) to find the frequency and
chirp mass distribution, which is beyond the scope of this paper. However, the
uncertainty in the orbital frequency of a binary on the detection threshold
(${\rm SNR}=7$) for a 4-year LISA mission is $2.5\times 10^{-9}\,\mathrm{Hz}$
and $1.0\times 10^{-9}\,\mathrm{Hz}$ for a 10-year mission (calculated using
Eq. C4). Therefore, we expect that the SKA could likely isolate the correct
binary pulsar to match to a LISA detection even when several are present in
the sky localisation region.

###  5.3 Caveats

Our predictions are subject to various uncertainties which can be broadly
divided into two different categories: those arising from the progenitor
models for the population of DCOs and those arising from the choices we have
made when placing these DCOs in our model for the Milky-Way. Although we are
unable, at present, to evaluate the impact of all these uncertainties, the
reader should nevertheless keep in mind that they are likely very substantial.
Most of these concerns are not unique to these study, but apply to most of the
predictions available in present literature. We highlight a few main concerns.

##### Progenitor models

Our binary-star progenitors models have been computed with a rapid population
synthesis code (see Sec. 2.1). This code relies on approximate parametric
prescriptions for the stellar evolutionary tracks of single tracks and simple
algorithms to mimic the effects of evolutionary and binary interaction
processes. Even though we explicitly consider the impact of some of the main
physics uncertainties (see Sect. A.3) the list of variations that we
considered is far from exhaustive. Moreover, it is by no means guaranteed that
the parametric prescriptions used in this code lead to realistic results, even
when varying the values of the parameters to their extremes. We stress in
particular the uncertainties affecting our most massive progenitor models.
Observational constraints are scarce for high mass stars and practically non-
existent for the rapid evolutionary phases (e.g. Langer, 2012; Mapelli, 2021).
This is even more true for the evolution of massive stars at low metallicity.
In addition to our limited understanding of massive stars, we note that the
rapid population synthesis code, such as the one employed to compute the
models used in this study, rely on extrapolations of the original fitting
formulae to approximate the evolutionary tracks for these higher mass
progenitors (Hurley et al., 2000, 2002).

##### Populating the Milky Way

Our Milky Way model is semi-empirical and has been calibrated based on
observations. Unfortunately, the early evolution of the (metallicity
dependence of the) star formation history is poorly constrained. We do not
expect this to be a major concern, as most of the double compact objects have
relatively short delay times of less than $2\,\mathrm{Gyr}$ (see Fig. 3e), but
this is a caveat that should be kept in mind. Furthermore, to estimate the
rate of detectable systems, we rely on normalisation choices (e.g. how many
detectable double compacts are formed per unit of star formation). This
depends heavily on the initial mass function, as low mass stars account for
most of the mass while high mass stars are the progenitors of double compact
objects. Further choices, such as the binary fraction and the initial
distributions of binary parameters also play a lesser but probably still
significant role de Mink & Belczynski (e.g. 2015); Chruslinska et al. (e.g.
2017); Klencki et al. (e.g. 2018).

We also note that, for reasons of computational efficiency, we have not
accounted for the spatial velocities resulting from the Blaauw-Boersma kick
(Blaauw, 1961; Boersma, 1961). In test simulations we find that accounting for
this spreads out the population (increasing the typical height above the
Galactic plane and Galactocentric radius), but we find that the impact on the
rate is limited. In light of the other much larger uncertainties, we felt that
this was justified (see however, e.g., Brandt & Podsiadlowski 1995; Abbott et
al. 2017b). We have further ignored a possible contribution coming from the
Galactic halo, as Lamberts et al. (2018) estimates this not be significant.
However, this may not be true for other formation channels other than those we
have considered here.

###  5.4 Other formation channels

In this paper we considered the formation of NS and BH binaries formed via
isolated binary evolution, through the classical CE channel, the stable mass
transfer channel and variations on these (see Fig. 15). We did not consider
further possible contributions from other formation channels, which may play a
role.

We highlight the possible role of dynamical formation in globular clusters.
Kremer et al. (2018) predict, for a nominal 4-year LISA mission, that 21
sources will have SNR $>7$, of which 7 are BHBHs, 0 are BHNSs and 1 is a NSNS
(see Table 1 Kremer et al., 2018). This is significantly lower than the rates
we predict for nearly every model variation. If true, this would mean that
formation through isolated binary formation will dominate the LISA detections.

Banerjee (2020) investigates formation of LISA detectable BHBHs in young
massive and open stellar clusters and estimates approximately 128 BHBHs with
SNR $>5$ in a 5-year LISA mission (see Table 1, Column 9 Banerjee, 2020).
Although this is similar to the number we predict for our fiducial model, we
note these authors adopt a threshold SNR required for a detection that is
lower and a mission length is slightly longer than what is typically assumed
(i.e. SNR $>7$ and 4 years, as we have also adopted in our work). We expect
that, after correcting for this and making a fair comparison, our fiducial
model predicts about twice as many detections.

The contribution of triples systems (e.g. Antonini et al., 2017), or even
higher order multiple systems (e.g. Vynatheya & Hamers, 2021) will likely also
be of interest, in particular for the formation of eccentric sources. We are,
however, not aware of specific predictions for the detection rates that we can
compare to directly.

##  6 Comparison with previous studies

Figure 12: A table comparing previous studies of a similar nature to this
work. In the detection rate section darker shades of blue indicate higher
detection rates. In other sections, darker shades of green indicate more
detailed models or more accurate choices. Entries in detection rates with two
values listed correspond to the different mission lengths (penultimate
column). The works listed in the table are Nelemans et al. (2001), Belczynski
et al. (2010), Liu & Zhang (2014), Lamberts et al. (2018), Sesana et al.
(2020), Lau et al. (2020), Breivik et al. (2020) and Shao & Li (2021).  .

In Fig. 12, we compare our results to similar previous studies that
investigate the population of stellar-mass BHBHs, BHNSs and NSNSs that are
detectable with LISA. Fig. 12 details the expected detection rates predicted
by each paper as well as their assumptions regarding their Milky Way galaxy
model, binary population synthesis simulation and LISA mission specifications.
We only include papers that are similar to our work, such that they use
population synthesis and simulate sources in the Galactic plane. Moreover,
Fig. 12 does not include the numerous papers on the LISA WDWD population as we
do not make predictions for WDWDs.

##### Nelemans et al. (2001)

were the first to investigate the population of LISA detectable stellar-mass
double compact objects. We find a significantly higher detection rate for
BHBHs and BHNSs, as well as a slightly lower rate for NSNSs. We can understand
this difference from changes both to the specifications of LISA (such as the
mission length and SNR threshold for detection) and our understanding of
massive star evolution since the publication of their paper, which both
strongly affect the expected detections rates.

##### Belczynski et al. (2010)

built upon the work of Nelemans et al. (2001), by using a different population
synthesis code with two model variations and a multi-component model for the
Milky Way. They find a much lower detection rate for BHNSs and NSNSs (and
agreed on zero BHBHs) when compared to Nelemans et al. (2001). They state that
this discrepancy from Nelemans et al. (2001) comes from differences in their
population synthesis and an overall lower formation rate rather than any
changes to LISA detectability. The low total detection rate for all DCOs in
this paper compared to our work is unsurprising given the relatively high SNR
threshold of 10 and short mission length of 1 year. The reduced mission length
means that the source signal has much less time to accumulate, whilst also
fewer WDWDs can be resolved in this time, leading to a weaker signal and an
increased Galactic confusion noise relative to our work.

##### Liu & Zhang (2014)

performed a similar investigation using a different population synthesis code
and find higher rates than earlier works. Their lower detection threshold and
longer mission length compared to Belczynski et al. (2010) likely explains the
relatively increased rates. Yet their rates are still significantly below what
we find. This could be for several reasons; they assume all binaries are
circular both in their evolution and for detection. This means that systems
may not have inspiralled as far before the LISA mission or may appear to have
weaker gravitational waves when eccentricity is not accounted for. They also
use a simplified model for the Milky Way with a single disc of one metallicity
and constant star formation, whilst also using a mission length half what we
assume. Each of these factors likely contributes to the lower overall
detection rates.

##### Lamberts et al. (2018)

presented a new approach to the problem by using the FIRE simulation (Hopkins
et al., 2014) to distribute their sources rather than an analytical model of
the Milky Way and were the first paper in this area to incorporate metallicity
dependence into their Milky Way model. Sesana et al. (2020) followed up on
this paper using the same simulated BHBH population and presented updated
results for the number of expected BHBH detections. They find significantly
fewer BHBHs than our fiducial model despite using the same SNR threshold and
LISA mission length. The discrepancy between the results of Sesana et al.
(2020) and those presented in this work could be caused by different
treatments of eccentricity. Unlike our work, Sesana et al. (2020) assume that
all binaries are circular for the purpose of detection in LISA, which could
result in a lower number of detections by missing eccentric binaries that
appear as weaker signals when assumed to be circular. This is especially
relevant as we find that around $87\%{}$ of LISA detectable BHBHs are not
circular and around $21\%{}$ have significant eccentricity (see Fig. 3d). We
also improve upon this work by using a larger number of metallicity bins
compared to Sesana et al. (2020), since a low number of metallicity bins can
produce artificial features in the mass distribution of DCOs and possibly
affect the detection rate (see Dominik et al., 2015; Neijssel et al., 2019;
Kummer, 2020, and also appendix D for further discussion). Finally, it could
be that different implicit assumptions in their population synthesis code lead
to differences in our results (Toonen et al., 2014).

##### Lau et al. (2020)

focussed on the number of Galactic NSNS binaries that could be detected by
LISA. Their study uses the same population synthesis code, COMPAS, as this
work, though an earlier version. Despite this, their study finds a much larger
number of detections. They make several different physical assumptions in
their population synthesis, using the Fryer et al. (2012) rapid remnant mass
prescription, assuming the optimistic CE scenario, limiting the maximum
neutron star mass to $2\,\mathrm{M_{\odot}}$ and not implementing PISN.
However, we note that none of these assumptions strongly affect the NSNS LISA
detection rate (see bottom panel of Fig. 8, models K, L, M and O) and so this
is unlikely to entirely account for the differences. It is also important to
highlight that COMPAS has received several improvements and bug fixes since
Vigna-Gómez et al. (2018) (which contains the simulations used by Lau et al.
2020) and these could possibly have affected the formation rate of NSNSs. Yet
it is most likely that the remaining difference between our results is due to
way in which we simulate the Milky Way. Lau et al. (2020) use a model for the
Milky Way similar to that of Breivik et al. (2020), which we use to estimate
the impact of the choice of MW model in Appendix D. The Milky Way model by
Breivik et al. (2020) applies only two metallicity bins, while we consider a
range of metallicities between $10^{-4}$ and $0.03$. When applying the simpler
model for the Milky Way, we find that the NSNS detection rate is increased by
at least a factor of two. Lau et al. (2020) only uses a single metallicity,
and so assuming that all star formation happens at a single high metallicity
(which has a high efficiency of producing NSNS), could lead to an even greater
overestimate of the detection rate. Hence, we expect the low number of
metallicity bins in their Milky Way model to be the main driver behind the
discrepancy between our results.

##### Breivik et al. (2020)

introduced the population synthesis code COSMIC and presented detections for
many different DCO types in LISA using this code. They find that LISA will
detect 93 BHBHs, 33 BHNSs and 8 NSNSs in the Milky Way over a 4 year mission.
Breivik et al. (2020) make many physical assumptions that differ from our
fiducial model, the most notable being that they assume the optimistic CE
scenario and that case BB mass transfer is always unstable, whilst also using
a simpler model for the Milky Way (see Appendix D). Thus for better comparison
we ran our simulation using model F and the Milky Way model from Breivik et
al. (2020). This results in 97, 101 and 43 detections for BHBHs, BHNSs and
NSNSs respectively. Therefore, though we are in very good agreement for BHBHs,
we predict much higher rates for BHNSs and NSNSs. These differences are likely
due to using a different population synthesis code (COSMIC), which has
different underlying physics assumptions from COMPAS. Given our strong
agreement for BHBHs, it is possible that COSMIC and COMPAS handle NSs
differently and so lead to different detection rates for DCOs containing NSs.
However checking this would require a more in-depth study of the intrinsic
formation rate of DCOs containing NSs in the two codes.

##### Shao & Li (2021)

most recently investigated the detectability binaries containing BHs in LISA
using BSE and a relatively simple model for the Milky Way (assuming a uniform
flat disc, constant star formation and a single metallicity). They assume that
kicks for NSs formed through ECSN are slightly higher than our work
($50\,\mathrm{km}\,\mathrm{s^{-1}}$ instead of
$30\,\mathrm{km}\,\mathrm{s^{-1}}$). This may account for their particularly
low BHNS rate (as the binaries would be more likely to disrupt), which is a
factor of 20 lower than ours, but we expect their assumption of the optimistic
CE scenario, reduced Wolf-Rayet winds and lower SNR detection threshold could
partially offset this. As we show in Appendix D, their use of a simpler Milky
Way model, especially with only a single metallicity, would lead to an
underestimate of the BHBH and BHNS rates, which may explain the discrepancy in
our results.

Overall, since the work of Nelemans et al. (2001), in addition to the LISA
mission specifications, the methods that we use to simulate binaries and the
Milky Way have all changed significantly. We now predict that LISA detections
of these massive DCOs are dominated by BHBHs, rather than NSNSs, whilst the
absolute detection rates for BHBHs and BHNSs are much higher. Further studies
in this area could improve on this work by including the effects of systemic
kicks on the position of systems in the Milky Way and accounting for
contributions from other formation channels.

##  7 Conclusion & Summary

We provide predictions for the detection rate and population properties of
LISA detectable BHBH, BHNS and NSNS. We use a novel empirically-informed
analytical model for the metallicity dependent star formation history of the
Milky Way, calibrated against the APOGEE stellar spectroscopic survey. We use
this to model Monte-Carlo realisations of the present-day BHBH, BHNS and NSNS
populations in our Milky Way. For the binary population, we use the results of
a large grid of simulations performed with the rapid population synthesis code
COMPAS (Broekgaarden et al. 2021, Broekgaarden et al. (in prep.)). These
simulations have been optimised with the adaptive sampling algorithm
STROOPWAFEL to preferentially sample NS and BH binaries. In total these
comprise over two billion massive binaries that span 20 physics variations,
which represent the most common uncertainties in binary physics. To determine
the detectability of sources with LISA we use the LEGWORK package (Wagg et
al., 2021), that was specifically developed for this purpose and is publicly
available. We investigate the results expected for a 4- and 10-year LISA
mission. Our main conclusions can be summarised as:

  1. 1.

Total detections: We predict 30-370 detections in a 4-year LISA mission,
across all our simulations for varying physics assumptions. This increases to
about 50-550 for a 10-year LISA mission. Although the number of detections per
type can vary by about 2 orders of magnitude, we find that the total detection
rate is fairly robust, among the variations we have considered (see Table 3).

  2. 2.

Detections by type: For our fiducial model, we predict a total of $124\pm 11$
detections and out of these we find about $74{}\pm 9$ BHBHs, $42{}\pm 6$ BHNSs
and $8{}\pm 3$ NSNSs. The errors quoted here are the $1$-$\sigma$ Poisson
uncertainties resulting from the random initialisation of the Milky Way (see
Table 3).

  3. 3.

Physics variations: Among the model variations we consider, we find that our
predictions for the rates for the different DCO types are robust within a
factor of 2 of the fiducial rate, with the following exceptions. For BHBHs,
the rate is most sensitive to the treatment of common envelope phases or an
increase of the WR wind mass loss. For BHNSs and NSNSs, the assumptions
regarding the assumed common envelope ejection efficiency, treatment of case
BB mass transfer and the kicks are most important. In addition, the assumed
mass transfer efficiency impacts the BHNS (see Fig. 8).

  4. 4.

Probing the black hole mass distribution and the lower mass gap: We expect
LISA to predominantly detect lower mass BHs (with 90% of BHBH and BHNSs having
BH masses lower than $11\,\mathrm{M_{\odot}}$ in our fiducial simulations) in
stark contrast to current ground-based detectors which are heavily biased
towards high mass systems. For our fiducial simulation, we predict that
approximately 69 systems with a component with a mass between
$2.5$-$5\,\mathrm{M_{\odot}}$ would be detected by a 4-year LISA. This implies
that LISA can potentially make important contributions to the debate about the
existence of a lower mass gap (see Fig. 9).

  5. 5.

Eccentricity distribution: We find that for all DCO types a large fraction of
detectable systems still have nonzero eccentricities ($e>0.01$) when entering
the LISA band, which can be used to distinguish them from the more numerous
WDWD binaries, which are largely expected to be circular. In particular, for
our fiducial model, we find that this is the case for around three quarters of
detectable binaries. Furthermore, around 16% of detectable binaries have
eccentricities that are so high ($e>0.3$) that the emission at frequencies
corresponding to higher order harmonics start to dominate (see Fig. 3).

  6. 6.

Distinguishing from WDWD sources: For about half of all detections we expect
that we will be able to confidently determine the type of compact objects
involved and this increases to 60% for a 10-year LISA mission (see Sec.
5.1.1).

  7. 7.

Chirp mass determinations: For about 10% of systems we expect to be able to
determine the chirp mass better than 10% and this increases to 15% for a
10-year LISA mission (see Fig. 6).

  8. 8.

Prospects for finding EM counterparts: We expect about 13% of detections with
a sky localisation better than 1 degree for our fiducial model (though the
fraction remains roughly constant among model variations). This fraction
remains the same for a 10-year LISA mission, meaning that the number increases
proportionally. This will be of interest for electromagnetic searches for
counterparts, in particular for radio pulsar searches with SKA (see Fig. 5.2).

  9. 9.

Benefits of extending the LISA mission: The number of detections scale
approximately as $\sqrt{T_{\rm obs}}$, where $T_{\rm obs}$ is the mission
length. Therefore, extending the LISA mission from 4- to 10-years increases
the number of detections by about 60% for each model variation. A further
important benefit is the improvement of the characterisation of the sources,
since the relative error on the frequency derivative (which dominates the
relative error in the chirp mass) scales as $T_{\rm obs}^{-2.5}$ for
stationary sources (Eq. C5). We find that the number of systems with chirp
masses that can be measured better than 10% increases by a factor of 2.4 for
each model variation. In addition, the number of systems with a sky
localisation better than one degree increases by a factor of 1.5. Overall, the
number of sources that can be unambiguously distinguished from WDWDs increases
by almost a factor of 2 (see Section 5.1.1).

We thank Katie Brevik, Will Farr, Rob Farmer, Valeriya Korol, Floris Kummer,
Eva Laplace, Mike Lau, Tyson Littenberg, Ilya Mandel, Mathieu Renzo, Alberto
Sesana, Katie Sharpe, Simon Stevenson, the CCA GW group and COMPAS
collaboration for insightful discussions. TW also thanks Terence Lovatt for
his advice and help on improving an earlier draft of this project. Portions of
this study were inspired by the LISA Sprint at the Center for Computational
Astrophysics of the Flatiron Institute, supported by the Simons Foundation.
This project was funded in part by the National Science Foundation under Grant
No. (NSF grant number 2009131), the European Union’s Horizon 2020 research and
innovation program from the European Research Council (ERC, Grant agreement
No. 715063), and by the Netherlands Organization for Scientific Research (NWO)
as part of the Vidi research program BinWaves with project number 639.042.728.
We further acknowledge the Black Hole Initiative funded by a generous
contribution of the John Templeton Foundation and the Gordon and Betty Moore
Foundation.

## References

  * Abadie et al. (2010) Abadie, J., Abbott, B. P., Abbott, R., et al. 2010, Classical and Quantum Gravity, 27, 173001, doi: 10.1088/0264-9381/27/17/173001
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review X, 6, 041015, doi: 10.1103/PhysRevX.6.041015
  * Abbott et al. (2019) —. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2020a) —. 2020a, ApJ, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2017a) —. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, ApJ, 850, L40, doi: 10.3847/2041-8213/aa93fc
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020b, arXiv e-prints, arXiv:2010.14527.  https://arxiv.org/abs/2010.14527
  * Abbott et al. (2020c) —. 2020c, ApJ, 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2021) —. 2021, ApJ, 913, L7, doi: 10.3847/2041-8213/abe949
  * Abdul-Masih et al. (2020) Abdul-Masih, M., Banyard, G., Bodensteiner, J., et al. 2020, Nature, 580, E11, doi: 10.1038/s41586-020-2216-x
  * Abt (1983) Abt, H. A. 1983, ARA&A, 21, 343, doi: 10.1146/annurev.aa.21.090183.002015
  * 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
  * Antoniadis et al. (2016) Antoniadis, J., Tauris, T. M., Ozel, F., et al. 2016, arXiv e-prints, arXiv:1605.01665.  https://arxiv.org/abs/1605.01665
  * Antonini & Rasio (2016) Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187, doi: 10.3847/0004-637X/831/2/187
  * Antonini et al. (2017) Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77, doi: 10.3847/1538-4357/aa6f5e
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * 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
  * Banerjee (2020) Banerjee, S. 2020, Phys. Rev. D, 102, 103002, doi: 10.1103/PhysRevD.102.103002
  * Barack & Cutler (2004) Barack, L., & Cutler, C. 2004, Phys. Rev. D, 69, 082005, doi: 10.1103/PhysRevD.69.082005
  * Bates et al. (2014) Bates, S. D., Lorimer, D. R., Rane, A., & Swiggum, J. 2014, MNRAS, 439, 2893, doi: 10.1093/mnras/stu157
  * Begelman et al. (1980) Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307, doi: 10.1038/287307a0
  * Belczynski et al. (2010) Belczynski, K., Benacquista, M., & Bulik, T. 2010, ApJ, 725, 816, doi: 10.1088/0004-637X/725/1/816
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016, Nature, 534, 512, doi: 10.1038/nature18322
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223, doi: 10.1086/521026
  * Belczynski et al. (2007) Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504, doi: 10.1086/513562
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Bertelli et al. (1994) Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&AS, 106, 275
  * Blaauw (1961) Blaauw, A. 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Boersma (1961) Boersma, J. 1961, Bull. Astron. Inst. Netherlands, 15, 291
  * Bolton (1972) Bolton, C. T. 1972, Nature, 235, 271, doi: 10.1038/235271b0
  * Bovy et al. (2019) Bovy, J., Leung, H. W., Hunt, J. A. S., et al. 2019, MNRAS, 490, 4740, doi: 10.1093/mnras/stz2891
  * Bovy et al. (2016) Bovy, J., Rix, H.-W., Schlafly, E. F., et al. 2016, ApJ, 823, 30, doi: 10.3847/0004-637X/823/1/30
  * Brandt & Podsiadlowski (1995) Brandt, N., & Podsiadlowski, P. 1995, MNRAS, 274, 461, doi: 10.1093/mnras/274.2.461
  * Breivik et al. (2020) Breivik, K., Coughlin, S., Zevin, M., et al. 2020, ApJ, 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Breivik et al. (2016) Breivik, K., Rodriguez, C. L., Larson, S. L., Kalogera, V., & Rasio, F. A. 2016, ApJ, 830, L18, doi: 10.3847/2041-8205/830/1/L18
  * Broekgaarden (2021a) Broekgaarden, F. S. 2021a, BHBH simulations from: Impact of Massive Binary Star and Cosmic Evolution on Gravitational Wave Observations II: Double Compact Object Mergers, 1, Zenodo, doi: 10.5281/zenodo.5651073
  * Broekgaarden (2021b) —. 2021b, BHNS simulations from: Impact of Massive Binary Star and Cosmic Evolution on Gravitational Wave Observations II: Double Compact Object Mergers, 1, Zenodo, doi: 10.5281/zenodo.5178777
  * Broekgaarden (2021c) —. 2021c, NSNS simulations from: Impact of Massive Binary Star and Cosmic Evolution on Gravitational Wave Observations II: Double Compact Object Mergers, 1, Zenodo, doi: 10.5281/zenodo.5189849
  * Broekgaarden et al. (2021) Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, arXiv e-prints, arXiv:2103.02608.  https://arxiv.org/abs/2103.02608
  * Broekgaarden et al. (2019) Broekgaarden, F. S., Justham, S., de Mink, S. E., et al. 2019, MNRAS, 490, 5228, doi: 10.1093/mnras/stz2558
  * Brown (1995) Brown, G. E. 1995, ApJ, 440, 270, doi: 10.1086/175268
  * Chandrasekhar (1931) Chandrasekhar, S. 1931, ApJ, 74, 81, doi: 10.1086/143324
  * Chattopadhyay et al. (2020) Chattopadhyay, D., Stevenson, S., Hurley, J. R., Bailes, M., & Broekgaarden, F. 2020, arXiv e-prints, arXiv:2011.13503.  https://arxiv.org/abs/2011.13503
  * Chruslinska et al. (2017) Chruslinska, M., Belczynski, K., Bulik, T., & Gladysz, W. 2017, Acta Astron., 67, 37
  * Collette et al. (2021) Collette, A., Kluyver, T., Caswell, T. A., et al. 2021, h5py/h5py: 3.3.0, 3.3.0, Zenodo, doi: 10.5281/zenodo.594310
  * Colpi et al. (2019) Colpi, M., Holley-Bockelmann, K., Bogdanovic, T., et al. 2019, arXiv e-prints, arXiv:1903.06867.  https://arxiv.org/abs/1903.06867
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61, doi: 10.1051/0004-6361/201527130
  * de Kool (1990) de Kool, M. 1990, ApJ, 358, 189, doi: 10.1086/168974
  * De Marco et al. (2011) De Marco, O., Passy, J.-C., Moe, M., et al. 2011, MNRAS, 411, 2277, doi: 10.1111/j.1365-2966.2010.17891.x
  * de Mink & Belczynski (2015) de Mink, S. E., & Belczynski, K. 2015, ApJ, 814, 58, doi: 10.1088/0004-637X/814/1/58
  * de Mink et al. (2009) de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243, doi: 10.1051/0004-6361/200811439
  * de Mink et al. (2007) de Mink, S. E., Pols, O. R., & Hilditch, R. W. 2007, A&A, 467, 1181
  * Dewdney et al. (2009) Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proceedings, 97, 1482, doi: 10.1109/JPROC.2009.2021005
  * Dewi et al. (2002) Dewi, J. D. M., Pols, O. R., Savonije, G. J., & van den Heuvel, E. P. J. 2002, MNRAS, 331, 1027, doi: 10.1046/j.1365-8711.2002.05257.x
  * Di Carlo et al. (2020) Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020, MNRAS, 498, 495, doi: 10.1093/mnras/staa2286
  * Dominik et al. (2015) Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, ApJ, 806, 263, doi: 10.1088/0004-637X/806/2/263
  * Dosopoulou & Kalogera (2016a) Dosopoulou, F., & Kalogera, V. 2016a, ApJ, 825, 70, doi: 10.3847/0004-637X/825/1/70
  * Dosopoulou & Kalogera (2016b) —. 2016b, ApJ, 825, 71, doi: 10.3847/0004-637X/825/1/71
  * du Buisson et al. (2020) du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941, doi: 10.1093/mnras/staa3225
  * Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126, doi: 10.1038/340126a0
  * El-Badry & Quataert (2020) El-Badry, K., & Quataert, E. 2020, MNRAS, 493, L22, doi: 10.1093/mnrasl/slaa004
  * Eldridge et al. (2020) Eldridge, J. J., Stanway, E. R., Breivik, K., et al. 2020, MNRAS, 495, 2786, doi: 10.1093/mnras/staa1324
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103, doi: 10.1088/0004-637X/741/2/103
  * Finn & Thorne (2000) Finn, L. S., & Thorne, K. S. 2000, Phys. Rev. D, 62, 124021, doi: 10.1103/PhysRevD.62.124021
  * Frankel et al. (2018) Frankel, N., Rix, H.-W., Ting, Y.-S., Ness, M., & Hogg, D. W. 2018, ApJ, 865, 96, doi: 10.3847/1538-4357/aadba5