# Predicting gravitational wave signals from BPASS White Dwarf Binary and Black Hole Binary populations of a Milky Way-like galaxy model for LISA

Petra Tang, Jan Eldridge, Renate Meyer, Astrid Lamberts, Guillaume

Boileau, Wouter van Zeist

► To cite this version:

Petra Tang, Jan Eldridge, Renate Meyer, Astrid Lamberts, Guillaume Boileau, et al. Predicting gravitational wave signals from BPASS White Dwarf Binary and Black Hole Binary populations of a Milky Way-like galaxy model for LISA. Monthly Notices of the Royal Astronomical Society, 2024, 534 (3), pp.1707-1728. 10.1093/mnras/stae2154. hal-04609744

HAL Id: hal-04609744

https://hal.science/hal-04609744v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Predicting gravitational wave signals from BPASS white dwarf binary and black hole binary populations of a Milky Way-like galaxy model for LISA

P. Tang<sup>1,1*</sup> J. J. Eldridge<sup>1,1*</sup> R. Meyer,<sup>2,*</sup> A. Lamberts<sup>1,3,4</sup> G. Boileau<sup>1,5</sup> and W. G. J. van Zeist<sup>1,6,7</sup>

Accepted 2024 September 13. Received 2024 September 12; in original form 2023 September 25

## ABSTRACT

Galactic white dwarf binaries (WDBs) and black hole binaries (BHBs) will be gravitational wave (GW) sources for LISA. Their detection will provide insights into binary evolution and the evolution of our Galaxy through cosmic history. Here, we make predictions of the expected WDB and BHB population within our Galaxy. We combine predictions of the compact remnant binary populations expected by stellar evolution from the detailed Binary Population and Spectral Synthesis (BPASS) code, with a Milky Way analogue galaxy model from the Feedback in Realistic Environment (FIRE) simulations. We use PHENOMA and LEGWORK to simulate LISA observations. Both packages make similar predictions that on average four Galactic BHBs and 673 Galactic WDBs are above the signal-to-noise ratio (SNR) threshold of 7 after a four-year mission. We compare these predictions to earlier results using the binary star evolution (BSE) code with the same FIRE model galaxy. We find that BPASS predicts a few more LISA observable Galactic BHBs and a twentieth of the Galactic WDBs. The differences are due to the different physical assumptions that have gone into the binary evolution calculations. These results indicate that the expected population of compact binaries that LISA will detect depends very sensitively on the binary population synthesis models used and thus observations of the LISA population will provide tight constraints on our modelling of binary stars. Finally, from our synthetic populations, we have created mock LISA signals that can be used to test and refine data processing methods of the eventual LISA observations.

**Key words:** Galaxy: stellar content – gravitational waves – (stars:) white dwarfs – stars: black holes – binaries: general.

## 1 INTRODUCTION

The study of gravitational wave (GW) signals in the mHz band will help researchers reveal characteristics of the processes of binary evolution and the history of our Galaxy. The observations of the future Laser Interferometer Space Antenna (LISA) mission will provide the first results in this observational window and will set the stage for a new era of astrophysics. Led by the European Space Agency (ESA), LISA is designed to detect GW signals in the frequency ranges from $10^{-4}$ to 1 Hz. LISA is still in the development phase and is currently scheduled to launch in 2035. While the scientific community eagerly awaits its deployment, we have turned to simulations and predictive modelling to gain a glimpse of the signals LISA will eventually observe. By simulating data and predicting the GW signals that LISA is expected to detect, we can lay the groundwork for future discoveries and advance our

understanding of the cosmos. In particular, LISA will be able to detect Galactic binaries and resolve many individual ultracompact binaries (UCBs, e.g. Amaro-Seoane et al. 2023). These resolved UCBs mostly consist of white dwarf binaries (WDBs, e.g. Nelemans, Yungelson & Portegies Zwart 2004; Nissanke et al. 2012; Korol et al. 2017, 2022; Kremer et al. 2017; Lamberts et al. 2019; Long et al. 2020; Breivik et al. 2020a) and a smaller number of neutron stars and/or stellar black hole binaries (BHBs; e.g. Hils, Bender & Webbink 1990; Nelemans, Yungelson & Portegies Zwart 2001b, 2004; Lamberts et al. 2018; Andrews et al. 2020; Lau et al. 2020; Wagg et al. 2022d).

Galactic WDBs have been recognized as the primary contributor to signals detected by space-based GW observatories (Hils et al. 1990). However, there are three main source classes, namely the UCBs, massive binary black holes (MBHBs), and extreme mass ratio inspirals (EMRIs) as described within Amaro-Seoane et al. (2023). The majority of UCB sources are likely too distant to be individually detected by LISA. Instead, they will contribute a significant stochastic GW signal in LISA's frequency band originating from the millions of unresolved WDBs systems in the Galaxy (Farmer & Phinney 2003; Rosswog et al. 2009). This is sometimes referred to as confusion

<footnote n="1"><sup>1</sup>Department of Physics, University of Auckland, Private Bag 92019, Auckland, New Zealand</footnote>

<footnote n="2"><sup>2</sup>Department of Statistics, University of Auckland, Private Bag 92019, Auckland, New Zealand</footnote>

<footnote n="3"><sup>3</sup>Observatoire de la Côte d'Azur, Université Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, F-06304 Nice, France</footnote>

<footnote n="4"><sup>4</sup>Observatoire de la Côte d'Azur, Université Côte d'Azur, CNRS, Artemis, Bd de l'Observatoire, F-06304 Nice, France</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, University of Antwerpen, Prinsstraat 13, 2000 Antwerpen, Belgium</footnote>

<footnote n="6"><sup>6</sup>Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL, Nijmegen, the Netherlands</footnote>

<footnote n="7"><sup>7</sup>Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, the Netherlands</footnote>

<footnote>* E-mail: tangxiaocao@gmail.com (PT); j.eldridge@auckland.ac.nz (JJE); renate.meyer@auckland.ac.nz (RM)</footnote>
noise. However, resolvable sources will stand out as distinct signals above the instrumental and confusion noise. It is expected there will be approximately 10⁴ resolved WDBs, as predicted by many studies such as Nelemans et al. (2001a), Farmer & Phinney (2003), Ruiter et al. (2010), Korol et al. (2017), Lamberts et al. (2019), and Breivik et al. (2020a).

Although possible, detecting the stellar mass BHB population is a challenge; this is because they are rarer than WDBs and there are no known examples in our Galaxy as the only way to detect them will be with LISA. Until now, our knowledge comes from X-ray binaries and GW transients arising from the merger of BHBs (Casares, Jonker & Israeliian 2017; Abbott et al. 2021). The electromagnetically identified high-mass X-ray binaries that may eventually form BHBs in the Milky Way (MW) and nearby galaxies are also potential LISA sources (Belczynski, Bulik & Fryer 2012; Belczynski et al. 2013; van den Heuvel 2019). While the BHB that cause GW transients would have been detectable by LISA before merging, if they had resided within our Galaxy.

Comprehending the Galactic sources of mHz GWs is essential to both understand the sources themselves but also to be able to then remove them to detect possible cosmological sources and other populations such as EMRIs (Cornish et al. 2021). A binary population synthesis model is required to aid this search. In this work, we aim to use results from our own binary population synthesis results to predict the GW signals and understand what the signals are expected to be. At the same time, we contrast and compare our predictions to the predictions of Lamberts et al. (2018, 2019), who used the same galaxy model, but a different binary population synthesis code. Thus, we can explore how predictions vary between population synthesis codes, and explore the reasons for differences between these predictions.

The work of Lamberts et al. (2018, 2019), presented a method for predicting the population of BHBs and WDBs in the MW using cosmological simulations with the Feedback In Realistic Environments<fnref n="1" /> (FIRE) model (Hopkins et al. 2014). These FIRE simulations deal with the formation and evolution of stars and galaxies over the history of the Universe. The simulations provide information on how many stars form through the age of the Universe and how their metallicity varies. They attached a binary population synthesis model to the FIRE results, creating synthetic DWD and BHB populations, by matching stellar population predictions to the FIRE simulation star formation, ages, and metallicities.

In our study, we also use the same FIRE model to create a MW-like galaxy that is occupied by compact binaries, and then predict GW signals from these populations. Most binary population synthesis codes, such as the one that was implemented in Lamberts et al. (2018, 2019), are rapid population synthesis codes, which are known to be efficient. However, in our study, we use binary populations computed by the Binary Population and Spectral Synthesis (BPASS<fnref n="2" />) v2.2.1 code (Eldridge et al. 2017; Stanway & Eldridge 2018; Eldridge & Stanway 2022), which models the evolution of binary stars and stellar structure in detail rather than using analytic fits to single star models. This leads to different outcomes, especially increasing the stability of mass transfer (see Briel, Stevance & Eldridge 2023), which leads to differences in the predicted Galactic binary populations and the subsequent GW signals. By contrasting and comparing the predictions of the two sets synthetic GW signals from the BPASS

and BSE binary populations, we gain insights into the physics of interacting binaries.

This paper outlines the importance of detailed binary population synthesis codes for predicting GWs in LISA frequency band. In Section 2, we describe how we create our Galactic binary populations, and examine the details of our simulated WDB and BBH populations. In Section 3, we present the predicted GW signals for our synthetic populations and determine how many binary systems may be observed by LISA. In Section 4, we take the next step to calculate mock modulated GW signals from LISA for our different binary populations. Then in Section 5, we discuss our results, comparing to other predictions by other studies, before presenting our conclusions in Section 6.

## 2 FIRE AND BPASS

Our simulated predictions serve as crucial stepping stones, preparing researchers for the groundbreaking observations that LISA will provide once it is launched. This study combines a realistic model for the star formation history (SFH) and metallicity evolution for a Milky Way-like galaxy (Section 2.1), and a detailed binary population synthesis model, BPASS (Section 2.2).

### 2.1 FIRE simulation

To be able to compare our work to the studies by Lamberts et al. (2018, 2019), we use the same MW-like galaxies from the FIRE suite (Hopkins et al. 2014). We use the same m12i results, known as the 'Latte' simulation (Wetzel et al. 2023). This simulation is based on the improved 'FIRE-2' version of the code from Hopkins et al. (2018) and ran with the code GIZMO<fnref n="3" /> (Hopkins 2015). GIZMO solves the equations of hydrodynamics using the mesh-free Lagrangian Godunov 'MFM' method. The analysis of the simulations is conducted using the publicly available PYTHON package GIZMO ANALYSIS (Wetzel & Garrison-Kimmel 2020). For more detailed description of how FIRE was combined with the BSE binary population code, see Lamberts et al. (2018, 2019).

### 2.2 BPASS

BPASS is a suite of computer programs that models the evolution of interacting binary stars. It combines individual stellar models to create synthetic stellar populations and make varied predictions such as the spectral energy distributions of individual stars or of the total stellar population, the rates of electromagnetic and gravitational transients and much more (see Eldridge et al. 2017; Stanway & Eldridge 2018; Eldridge, Stanway & Tang 2019; Briel et al. 2022, and references therein). The stellar evolution models incorporate the physical processes such as mass transfer (MT), stellar winds, and supernova explosions. The main unique feature of BPASS (first described in Eldridge, Izzard & Tout 2008) is that the stellar models are evolved in detail using a custom version of the Cambridge STARS code (Eggleton 1971). This allows the reaction of the stellar structure to binary interactions to be evolved in detail. This typically leads to MT that is more stable in our models than is commonly assumed in other population synthesis codes (Briel et al. 2023). Another feature of BPASS is the predictions have been widely tested against a large range of observations. This list is growing all the time, for example see Kobayashi et al. (2023) who study the impact of binary

<footnote n="1">¹https://fire.northwestern.edu/</footnote>

<footnote n="2">²http://bpass.auckland.ac.nz or http://warwick.ac.uk/bpass</footnote>

<footnote n="3">³http://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html</footnote>
interactions on chemical enrichment of galaxies using BPASS and several other binary population synthesis code results. In this work, we aim to take steps to make BPASS predictions on future LISA observations of WDBs and BHBs in our Galaxy.

Many other stellar population synthesis codes are rapid binary population synthesis codes (e.g. Hurley, Tout & Pols 2002; Breivik et al. 2020a; Riley et al. 2022). The investigations by Lamberts et al. (2018, 2019) used the binary stellar evolution (BSE) code (Hurley et al. 2002) where the massive binary star treatment had been modified. BSE follows the evolution of binary stars with the evolution approximated using analytic fits to single star models from Hurley, Pols & Tout (2000) and the binary evolution described in Hurley et al. (2002). This approach allows for very computationally efficient exploration of a wide range of parameters. BSE models, however, miss important effects in the stellar structure and MT compared with detailed binary population synthesis codes (e.g. Gallegos-Garcia et al. 2021; Klencki et al. 2021; Marchant et al. 2021; Briel et al. 2023). Detailed structure evolution calculations result in stable MT (SMT) being more common than found by rapid codes. SMT occurs when matter is transferred from one star to its companion without the orbit dramatically shrinking and the binary entering a phase of common-envelope evolution. SMT significantly influences binary stars evolution and is important for the formation of many types of binary compact objects. However, SMT mechanisms are complex and uncertain, with factors like accretion rates, angular momentum transfer, and interaction between stellar envelopes playing important roles. Briel et al. (2023) and Ghodla & Eldridge (2023) discussed the effects of SMT in BPASS models, and postulated that SMT with super-Eddington accretion is responsible for high-mass BH population with masses 35–100 $M_\odot$ as observed by the LIGO–VIRGO–KAGRA collaboration.

In this work, we use the stellar models from BPASS v2.2.1. The population of WDBs and BHBs has been calculated using the GW population synthesis code TUI (v2.jje), a GW binary population synthesis programme which is described in Ghodla et al. (2022), Stevance et al. (2023), Briel et al. (2023) and van Zeist, Eldridge & Tang (2023). The predictions are made using the fiducial BPASS initial binary parameters and initial mass functions. The initial mass function is model 135–300 described in table 1 in Stanway & Eldridge (2018), which is the IMF from Kroupa, Tout & Gilmore (1993) with an upper mass limit of 300 $M_\odot$. We use the initial binary parameter distribution of Moe & Di Stefano (2017). The kick used in the population synthesis is the Hobbs neutron star kick, which is adapted to a momentum kick for BHs (Hobbs et al. 2005). TUI takes the final binary parameters after both stars have a stellar remnant. It then evolves the binary using the equations of Peters (1964) to model the impact of gravitational radiation removing energy from the orbit. Unless the systems merge, we stop integration of the orbit at 100 Gyrs. The resultant model population is binned by age, frequency, and remnant masses at a resolution of 0.1 dex, while binary eccentricity has bin sizes of 0.1.

To combine the binary population with the FIRE galaxy, we build on our work in van Zeist et al. (2023) where we used these same TUI models to model individual stellar clusters. TUI outputs the numbers of WDBs and BHBs expected, with different masses, periods, eccentricities at different ages of the stellar population, per $M_\odot$ of star formation at a number of initial metallicities. This can then be combined with the FIRE simulation. The FIRE simulation is a list of ‘star’ particles each with an initial mass, age, and metallicity. We assume that each particle is equivalent to a star cluster. Thus, the mass, age, and metallicity are used to identify a parameter space in the TUI outputs to then estimate how many binary compact

<table><caption>**Table 1.** Number of Galactic binary systems in the MW-like galaxy predicted by two models in the frequency range of $10^{-5}$ to 1 Hz. The full predicted population from BPASS of frequency ranging from $10^{-10}$ to 1 Hz is given in parentheses.</caption><thead><tr><th>Events</th><th>BSE</th><th>BPASS</th></tr></thead><tbody><tr><td>BHB</td><td>25 335</td><td>9298<br>(495 586)</td></tr><tr><td>WDB</td><td>35 832 029</td><td>58 522 007<br>(272 789 350)</td></tr></tbody></table>

objects each star particle contains. For each star particle, we simply multiply the BPASS-predicted WDB and BHB numbers per stellar mass (which can be a fraction) by the initial star particle mass. Each integer value of binaries is then kept, the remainder is carried over to the next star particle rather than being lost. It is important to note that we work through the list of star particle in age order. Therefore, carried over values are associated with systems with similar age, thus the distribution of the eventual binary systems within the synthetic catalogue will be similar. This method means that we do not underestimate the total binary population; however, it might mean that we miss rare types of binary systems. The fraction of stars carried over is typically small, of the order of 0.1 for WDB or of the order of $10^{-8}$ for BHB. Another way to evaluate the importance of the carry over is that there are 13 976 485 star particles in the galaxy model we use. This means for the full frequency range of binaries, on average for BPASS, each star particle hosts 19.52 WDBs and 0.35 BHBs. Given that the particles have the same age, spurious errors from carrying over the IMF weighting should be minor.

The final step is to randomly distribute the binary within the TUI output bin using a uniform random distribution. This avoids all the binaries having exactly the same parameters. This process results in a list of compact remnant binaries and their distribution within the synthetic galaxy. This allows for direct comparison to the predictions of Lamberts et al. (2018, 2019) who used the same synthetic galaxy but a different binary population synthesis code and different method of distributing star masses; the different ways of distributing star masses may not result in significant differences in mass distributions in the galaxy; however, this needs to be verified. The aim of our study is to provide insight into the importance of how the binaries are modelled for predicting the LISA observable population.

## 2.3 Comparing the BSE and BPASS synthetic populations in FIRE galaxy m12i

Having matched the FIRE m12i star particles with the binary population synthesized by BPASS, we can now compare to the BSE models of Lamberts et al. (2018, 2019). Both sets of calculated outputs include Cartesian galactic coordinates, orbital frequency, primary and secondary masses, age and metallicity of each binary. These outputs are used to create synthetic LISA observations in Sections 3 and 4.

### 2.3.1 The number of Galactic WDB and BHB

The numbers of Galactic binary systems predicted by the two binary population synthesis models are shown in Table 1. There are fewer BHB and more WDB systems predicted by BPASS than Lamberts et al. (2018, 2019), respectively. However, the number of WDB systems in both samples is closer in magnitude than the BHB populations. Most of these binaries are expected to be less than
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Large-scale 2D number density of BHB populations in the simulated galaxy from BPASS and BSE model populations viewed face on. These are the BHBs that are within the with frequencies > 10⁻⁵ Hz. BHBs are present in the galactic bulge, discs, galactic halo, and a satellite galaxy in both model populations. The population density is also highest at the centre of the galaxy for both models.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Large-scale 2D number density of WDB populations in the simulated galaxy from BPASS and BSE model predictions viewed face on. WDBs are present in the galactic bulge, discs, galactic halo, tidal streams, and satellite galaxies in both model predictions. The population density is highest at the centre of both models.</figcaption></figure>

20 kpc from the Galactic centre (cf. Seto 2002; Nelemans et al. 2004; Lamberts et al. 2019). However, because of the Galactic halo and tidal streams, some objects are at greater distances of up to 300 kpc, as shown in Figs 1 and 2.

### 2.3.2 Spatial density of Galactic WDBs

We have attempted to verify if our WDB model reproduces the expected space density of WDBs around the Sun. In Fig. 3, we show the expected number of WDBs and BHBs as we increase the distance of the objects included with the observed volume. The dotted line represents the Double Degenerate WDs (DDWDs) space density prediction, that are extrapolated from the DDWD population within 25 pc of the Sun surveyed by Holberg et al. (2016). Here, we use these DDWDs to verify the space density of the modelled WDBs. The space density of the WD population is estimated to be 4.8 ± 0.5 × 10⁻³ pc⁻³ at the Sun’s location. Within this WD population, 26 per cent are binaries, and 10 per cent are DDWDs. Using this estimated observed space density, we can extrapolate how many WDBs should be within the Galaxy. In Fig. 3, we can see that the BPASS WDB predictions are closer to the number of DDWDs from the extrapolation of the local space density. While BSE is an order of magnitude below the extrapolation. Our BPASS predictions includes the full frequency range of orbital periods (10⁻¹⁰–1 Hz), this is the main reason in the difference between BPASS and BSE and why the former matches the extrapolated DDWD space density within 10 kpc. We note that this extrapolation of the observed space density is approximate and does not take into account the Galactic disc structure, and the turnover in the theoretical numbers beyond 10 kpc demonstrates this within Fig. 3.

### 2.3.3 The total binary mass distribution

A summary of the distribution of total binary masses for the synthetic binary populations is shown in Table 2. The total mass of the BSE BHB population has higher median and mean values compared with the BPASS BHB population; and the maximum BHB in BPASS is more than twice the maximum BHB mass from BSE models. The total binary mass of BPASS WDB population has a lower mean and higher median values than those of the BSE; more low-mass WDBs survive in BPASS than in BSE.

The BPASS and BSE codes have differences in their assumed binary evolution models. This leads to different evolutionary outcomes for stars of the same initial mass, which here can be seen
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Cumulative number distribution of binary systems at different distances from the Sun. Blue and coloured lines represent Galactic BHB and WDB populations with frequencies above $10^{-5}$ Hz, respectively. The BPASS models are presented by solid lines and the BSE models by dashed lines. The solid green line represents BPASS WDBs with full frequency range. The observationally driven extrapolated DDWDBs population prediction is represented by the dotted straight line.</figcaption></figure>

<table><caption>**Table 2.** Summary of the mean, median, minimum, maximum, and standard deviation values of the total binary mass distribution of each WDB and BHB population with frequencies between $10^{-5}$ and 1 Hz. The full predicted population from BPASS of frequency ranging from $10^{-10}$ to 1 Hz is given in parentheses.</caption><thead><tr><th rowspan="2"></th><th colspan="2">BHB</th><th colspan="2">WDB</th></tr><tr><th>BSE</th><th>BPASS</th><th>BSE</th><th>BPASS</th></tr></thead><tbody><tr><td>Maximum</td><td>85.40 M<sub>☉</sub></td><td>191.10 M<sub>☉</sub><br>(249.60 M<sub>☉</sub>)</td><td>2.39 M<sub>☉</sub></td><td>2.38 M<sub>☉</sub><br>(2.82 M<sub>☉</sub>)</td></tr><tr><td>Median</td><td>25.53 M<sub>☉</sub></td><td>19.91 M<sub>☉</sub><br>(25.98 M<sub>☉</sub>)</td><td>0.90 M<sub>☉</sub></td><td>1.03 M<sub>☉</sub><br>(1.02 M<sub>☉</sub>)</td></tr><tr><td>Std. dev.</td><td>13.3 M<sub>☉</sub></td><td>19.4 M<sub>☉</sub><br>(20.81 M<sub>☉</sub>)</td><td>0.34 M<sub>☉</sub></td><td>0.29 M<sub>☉</sub><br>(0.27 M<sub>☉</sub>)</td></tr><tr><td>Mean</td><td>28.17 M<sub>☉</sub></td><td>26.21 M<sub>☉</sub><br>(30.10 M<sub>☉</sub>)</td><td>1.03 M<sub>☉</sub></td><td>1.07 M<sub>☉</sub><br>(1.07 M<sub>☉</sub>)</td></tr><tr><td>Minimum</td><td>6.72 M<sub>☉</sub></td><td>5.18 M<sub>☉</sub><br>(4.50 M<sub>☉</sub>)</td><td>0.39 M<sub>☉</sub></td><td>0.32 M<sub>☉</sub><br>(0.22 M<sub>☉</sub>)</td></tr></tbody></table>

in the predicted masses of the WDB systems. In this, we consider the primary star to always be the initially more massive star. MT in a binary occurs at different points of evolution, determined by the binary physics included as well as the initial masses and separations of the stars. While some of the differences are shown in Table 2, they are also apparent in the distribution of the masses of each remnant in a binary. We show in Figs 4 and 5, the distribution of the masses for both remnants in the of BHBs and WDBs, respectively. For the BHBs, we see a wider range of BHB masses, as well as examples of the most massive BH coming from the initially less massive star in BPASS results. This is also apparent for the WDBs where for BPASS the majority of the most massive WDs are from the initially less mass star.

For the BHBs, we expect the broader range of BH masses for BPASS is inherent in the stellar evolution models. If we compare the initial-to-final mass relations from the codes, using fig. 2 of

Lamberts et al. (2018) and fig. 17 of Eldridge et al. (2017), the range of black hole masses from a certain initial mass is wider for BPASS than BSE. This is likely due to the detailed stellar-evolution models that are calculated for BPASS populations. The mass-loss in binary interactions has a more significant impact on core masses and thus remnant masses than found in rapid codes. The impact of binary interactions on the remnant masses in detailed stellar evolution models has also been investigated by Laplace et al. (2021) and Patton, Sukhbold & Eldridge (2022).

For the BHBs, this mass reversal is relatively rarer than for WDBs. We find that this is because super-Eddington accretion onto the first BH is formed, it remains the most massive BH for most binaries as discussed by Briel et al. (2023). While for the WDBs, we find in BPASS that MT for low-mass stars is more stable than for BSE with more mass exchange from the primary to secondary star with fewer occurrences of common-envelope evolution (CEE). Especially during the first occurrence of Roche lobe Overflow (RLO) before the primary star becomes a WD. The reason for this stability is that BPASS uses detailed stellar evolution models for the donor star so MT stability is determined by how the star's radius responds to mass-loss. If the mass-loss stops, the radial expansion of the donor than mass transfer is stable. While if it continues to grow then a CE phase occurs. For BSE, as it is a rapid code, limits on stability are assumed determined by the evolutionary state of star and the mass ratio at onset of MT (e.g. Hurley et al. 2002). This may overestimate the number of systems that experience unstable MT. Such increased MT stability has been found by others who have used detailed stellar evolution models to follow binary interactions (e.g Gallegos-Garcia et al. 2021; Klencki et al. 2021; Marchant et al. 2021; Briel et al. 2023). This higher occurrence of MT stability is the leading cause of the differences in predicted GW signals for WDBs between BPASS and BSE. We find that 93 per cent of BPASS WDBs in the LISA band have the most massive WD from the secondary star, compared to only 1.2 per cent for the BSE models. The difference is less extreme for the BHBs with only 7.0 per cent of BPASS secondary stars being the most massive BHs, with 4.5 per cent for the BSE model.

### 2.3.4 Galactic map of Galactic WDBs and BHBs

Our next step is to understand whether LISA will be able to observe these WDBs and BHBs. To this end, we assume that the LISA constellation's position is at the Sun's location, at a distance of 8 kpc from the model galaxy's centre. Thus, the distance to each of the sources is computed from this point. The large-scale density maps of the Galactic binary population densities are shown in Figs 1 and 2, where all colour bars represent numbers of binary systems in each bin. The inner 30 kpc for the Galactic centre are presented in appendix Figs D1 and D2.

On the large scale map of 100 kpc radius about the LISA constellation shown in Fig. 1, both galaxies contain a visible satellite galaxy. On a small-scale, shown in Fig. D1, around a 25 kpc radius about the Sun's position, there are noticeable differences between the two BHB galaxies. BPASS BHBs are found mostly near the bulge, and Lamberts et al. (2018) BHBs spread out more evenly from the bulge towards the edge of the Galactic disc and halo, this can be seen in Fig. 3 as well. The distributions over the in the inner 10 kpc of the Galaxy are similar but with BSE predicting many more BHBs than BPASS. However, beyond 10 kpc the difference between BSE and BPASS becomes greater. With BSE predicting many more BHBs beyond this distance, as can be seen in Fig. 3. This difference is because we begin to sample the galactic halo where there are more
<figure><img src="image_5.png" /><figcaption>Figure 4. The distribution of primary M₁ and secondary M₂ masses in BHBs from BPASS (left) and BSE (right). The linear line represents M₁ = M₂, where above the line M₁ < M₂ and below the line M₁ > M₂.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. The distribution of primary M₁ and secondary M₂ masses in WDBs from BPASS (left) and BSE (right). The linear line represents M₁ = M₂, where above the line M₁ < M₂ and below the line M₁ > M₂.</figcaption></figure>

low-metallicity stars so the difference between BSE and BPASS in this regime become apparent.

Fig. 2 (small-scale map in appendix D2) shows the distribution of WDB systems. Both galaxies appear to be similar in structure. Furthermore, both galaxies contain visible satellite galaxies and tidal streams. These satellite galaxies' sizes are similar to the Small Magellanic Cloud in the MW galaxy.

### 2.3.5 $M_c$ distribution

In GW observations, a key parameter for any binary is the chirp mass $M_c$. It is a combination of the component mass and is defined as

$$M_c = \frac{(M_1 M_2)^{3/5}}{(M_1 + M_2)^{1/5}}, \qquad (1)$$

where M₁ and M₂ are primary and secondary masses of the binary, respectively. The amplitude of GWs is related to the frequency and $M_c$ of the binary source. Higher $M_c$ generally corresponds to more massive BHs and WDs and stronger GW signals. We find that BPASS predictions for BHB masses are higher, shown in Table 2. While in Fig. 6, we can see that only BPASS results have BHB with $M_c$ above 40 M☉. One reason for this is related to the initial mass ranges of stars used. In Lamberts et al. (2018), the initial mass of the stars has a maximum of 120 M☉ with the assumption more massive stars all

experienced Pair-Instability Supernova (PISNe) and so produce no black holes. However, in BPASS, the maximum initial stellar mass used is 300 M☉ and PISNe are assumed to occur when the star ends its evolution with a helium-core mass in the range between 64 and 133 M☉ thus preventing black holes with high masses from forming. As we note above, binary interactions in detailed evolution models have a more significant impact on the core masses and thus eventual fate of stars. Thus, while BHs from stars above 120 M☉ are rare, they do occur. A second, and the more significant reason, is the super-Eddington accretion assumption within BPASS populates this upper range of masses as detailed in Briel et al. (2023) (see also van Son et al. 2020; Ghodla & Eldridge 2023). Thus, these differences together explain the significant difference in the maximum chirp masses of the BHB predictions.

Next, if we look at chirp masses below 40 M☉ shown in Fig. 6, the general shape of the BPASS and BSE distributions are similar, where lower mass binaries present in with higher frequencies. With the majority of BHBs being in a $M_c$ range from 3 to 10 M☉, the BSE BHB population is more numerous than the BPASS BHB population, and has multiple systems with low $M_c$ and low frequencies. There are also more BPASS BHB systems present at the higher end of the LISA frequency band where LISA is most sensitive (10⁻³ Hz).

As discussed above the stability of MT in interacting binary stars is a key driver of the WDB mass distribution, affecting both the masses
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distribution of BHBs over $M_c$ and frequency, BPASS (left) and BSE (right). The BPASS BHB population has more systems with higher $M_c$ than BSE.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distribution of WDBs over $M_c$ and frequency, BPASS (left) and BSE (right). There are clear differences in the distributions between the two models that likely result from the different treatment of mass transfer and modelling of stellar structure.</figcaption></figure>

and the period distributions which are important for predicting the LISA signal from these objects. We can see the impact in the different distributions for BPASS and BSE in Fig. 7. There is a large triangular gap in the BSE distribution from $\log(M_c/M_\odot) = 0.3$ to 0.2 at frequencies of $10^{-4.5}$ to $10^{-3}$ Hz. While in the BPASS results, this region is populated. This is the result of the more stable model of MT in BPASS leading to more extremely low-mass WDs allowing for the broader range of possible WDB chirp masses. We find that distribution of chirp masses varies with age of the WDBs (see Fig. D5), where lower mass chirp masses are more likely to be older WDBs. Overall, there are fewer WDBs within the LISA frequency range. This is another result of the BPASS MT stability, because more stable MT means fewer CEE systems with their shorter periods. In addition, the BPASS CEE prescription is relatively weaker than these frequently used in BSE codes (see discussion in Stevance et al. 2023). Both of these mean that while we have more WDBs (because of fewer mergers in CEE), we have fewer WDBs with orbits in the observable range.

Finally, we also note that one factor that also changes close binary periods is magnetic-wind braking (e.g Rappaport, Verbunt & Joss 1983) for systems involving a single WD. As BPASS does not include this effect, in could be another reason why BPASS has fewer short period WDBs compared to BSE.

### 2.3.6 Metallicity distribution

The metallicity of a binary system affects its evolution and GW emission in several ways. Lower metallicity stars have weaker stellar winds and so can form more massive remnants at the end of their lives. We see this in Fig. 8 for the BPASS results, the BHB have higher masses at lower metallicities. Also, with weaker winds, stars can grow to larger radii so that binary interactions can become more common. Conversely, this means that it is possible binary interactions are less likely at high metallicity as the envelope will be removed before it can interact. Most of the BHB systems are found below solar metallicity $Z_\odot = 0.02$. BHB systems with $M_c$ greater than
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The distribution of BPASS BHBs over $M_c$ and metallicity of the progenitor stars.</figcaption></figure>

50 $M_\odot$ only exist in a narrow-metallicity band of around one-tenth solar metallicity. These most massive BHB systems above 50 $M_\odot$ are the result of SMT during the evolution of both the primary and secondary, with super-Eddington accretion onto the BHs leading to the most massive BHs.

We note that below a metallicity of $Z = 0.004$, we allow secondary stars in a binary to experience quasi-chemically homogeneous evolution. This is an atypical evolutionary pathway that is the result of mass transfer onto massive stars making the rotate so rapidly that rotationally induced mixing means they are fully mixed stars during their main-sequence evolution (Yoon & Langer 2005; Cantiello et al. 2007). This is assumed for stars at low $Z$ ($Z \le 0.004$), that are initially more massive than 20 $M_\odot$ and they accrete more than 5 per cent of their initial mass from the primary star. These stars are assumed to rapidly rotate and become fully mixed during their main sequence evolution. As such, they never fill their Roche lobe and donate material to the primary BH so it remains at its birth mass (see Briel et al. 2023, for further details). The implication here is that for such systems SMT onto the primary BH does not occur. This is the reason behind the horizontal linear feature at a $\log(Z/Z_\odot) = 10^{-0.5}$. Up to this metallicity, QHE does not occur, but below this where QHE is assumed to occur, SMT onto BHs only occurs for a few systems.

Below this metallicity, we see a much more gradual increase in the masses of black holes as metallicity decreases. This is similar to what we expect occurs in BSE. While we do not have the equivalent data for BSE, fig. 2 in Lamberts et al. (2018) allows us to infer that BSE has a regular increase in the BH masses as metallicity decreases.

The correlation between the metallicity and mass of WDs has important implications for detecting and characterizing GWs from WDB as well, as it can affect whether they can be observed by GW detectors. For both BPASS and BSE, the distributions are relatively smooth. However, for BPASS, we see a sudden shift at $\log(Z/Z_\odot) \approx 0.7$. We see a different trend in Fig. 9, the $M_c$ of WDBs also correlates with their progenitor systems’ metallicity. Specifically, higher metallicity systems tend to produce less massive WDs, leading to smaller cores and lower $M_c$ in resulting binary systems. This trend is also observed, although more weakly, in the BSE WDB population.

# 3 SYNTHETIC LISA OBSERVATIONS OF WDBS AND BHBS

## 3.1 Simulating LISA observations with PHENOMA and LEGWORK

With the population of WDBs and BHBs GW sources in an MW galaxy analogue, we can estimate the signals from these objects, understand what LISA will eventually observe and which objects will be resolvable. To do this, we make use of two publicly available packages for evaluating whether an object emitting GWs is detectable by LISA or not. These packages are PHENOMA (Ajith et al. 2007; Robson, Cornish & Liu 2019) and the LISA Evolution and Gravitational Wave ORbit Kit (LEGWORK, Wagg, Breivik & de Mink 2022b).

We note that PHENOMA has been updated to more generic and accurate models, such as the latest PHENOMB model (Hannam 2014; Schmidt, Ohme & Hannam 2015), which calculates spin-precession effects. However, PHENOMA, is sufficient for calculating strains and estimating signal-to-noise ratios (SNRs) in our study. While LEGWORK is being increasingly widely used by the community. By using both analysis packages, we are also able to evaluate if using different packages has any impact on evaluating how many sources in a binary population can be detected.

These packages work by first calculating the strain and amplitude spectral density (ASD), of a binary system based on its mass, orbital

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The distribution of BPASS (left) and BSE (right) WDBs over $M_c$ and metallicity of the progenitor stars.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The ASD for frequency range from $10^{-5}$ to 0.1 Hz of the BHB populations as calculated by PHENOMA. The solid line is the approximated full sensitivity curve combining the instrument noise and the Galactic confusion noise for a four-year LISA mission lifetime. The individual BHBs from BPASS population and BHBs from the BSE population are represented by points. The model populations are observed from the Sun's location at (8, 0, 0). The stars indicate the location of observed verification WDB and the circles are verification binaries of other types.</figcaption></figure>

period, and distance from LISA (see Ajith et al. 2007, Wagg et al. 2022b or Appendix A for details). This is then compared to the LISA sensitivity curve which has two components. The first of these is the pure instrumental noise of the LISA detector, the second is the confusion noise from unresolved WDBs. The GW ADS from a binary source can then be compared to evaluate its detectability, usually by assuming a minimum SNR. The details in our analysis are handled by PHENOMA and LEGWORK and we detail the noise curve and SNR calculation in Appendix B. Both PHENOMA and LEGWORK also output figures showing how the population of sources compare to the sensitivity curve as well as the number of systems that can be detected.

## 3.2 The GW emission of the synthetic binary populations

Before considering our mock observations of our synthetic WDB and BHB populations, we should consider what actual observations exist that can help us evaluate the validity of our predictions. A significant amount of observational work has been undertaken to find LISA verification binaries (e.g. Kupfer et al. 2018; Amaro-Seoane et al. 2023; Finch et al. 2023; Kupfer et al. 2023).

Verification binaries consist of a WD or neutron star primary and a compact subdwarf/WD/neutron star secondary that have been discovered by extant electromagnetic surveys. Most verification binaries have ultrashort orbital periods, less than an hour, producing GWs with frequency between $10^{-3}$ and $10^{-2}$ Hz. Only a fraction of the BHB and WDB populations from both BPASS and BSE models fall in this frequency range. The strains of many of the verification binaries are accurately predicted and calculated as their distances have been measured by Gaia (Gaia Collaboration et al. 2018). The expected SNRs of these verification binaries are all $\ge 5$ for a four-year LISA mission time, with many SNR values above the widely used LISA SNR threshold of 7. Detailed information on current verification binaries can be found in Kupfer et al. (2018,

2023), Amaro-Seoane et al. (2023), Finch et al. (2023), and the ADS library.<fnref n="4" />

To visualize the ASD of individual Galactic binary sources, we plot our model BHB and WDB populations, the 43 known LISA verification binary sources' ASD (marked as green dots and stars) and the LISA ASD noise sensitivity curves with and without confusion noise in Figs 10 and 11. The ASD for the sources is calculated by using equations (20–23) in Robson et al. (2019) with the sensitivity curve ASD from equations B1 to B7. Many verification binaries on these figures are well above the sensitivity curves, and thus, LISA is expected to detect many of these individual verification binary sources (Stroeer & Vecchio 2006; Burdge et al. 2020; Pelisoli et al. 2021). However, several are closer to the sensitivity curve and will be more challenging to resolve.

First, we notice most of the BHB population is below the sensitivity curve and will not be seen by LISA. Secondly, for both BSE and BPASS populations there are only a handful of BHBs above the sensitivity curve. BPASS models predict more BHBs with frequencies above $10^{-4}$ Hz than BSE, which can also be seen in Fig. 6. It is only systems above this frequency that can be detected by LISA, while BSE predicts more BHBs overall as BPASS reaches higher frequencies it predicts more detectable BHBs. We note that LEGWORK predicts almost exactly the same results as shown in Fig. D6. With the difference frequency distributions of the BSE and BPASS populations clearly evident.

Fig. 11 shows the ASD of the WDB populations for each of the two models against approximations of sensitivity curve ASD as calculated by PHENOMA. We note that the PHENOMA was incorrectly assuming systems with frequencies above $10^{-2}$ Hz were evolving during the LISA observations and would merge. We adjusted this part of the code so that all WDBs were treated as stationary sources. We notice that, in contrast to the BHB populations, there are orders of magnitude more WDBs above the LISA sensitivity curve which can

<footnote n="4">⁴https://ui.adsabs.harvard.edu/user/libraries/XOVGbNMAThua37oSEdV8LA</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** ASD versus frequency for the WDB populations as calculated by PHENOMA. The solid line is the approximated full sensitivity curve combining the instrument noise and the Galactic confusion noise for a four-year LISA mission lifetime. The individual WDBs from BPASS population and the BSE population are represented as points. The model populations are observed from the Sun's location at (8, 0, 0). The stars indicate the location of observed verification WDB and the circles are verification binaries of other types.</figcaption></figure>

potentially be individually resolved. There is also a great population of sources below the LISA sensitivity curve, the source of the confusion noise. We notice that there are more high-frequency and high-ASD ($> 10^{-18} \text{Hz}^{-1/2}$) systems in BSE model than the BPASS model, even though BPASS has more systems overall. The reverse of what is found for BHBs. Again, we have used LEGWORK to calculate the ASD as shown in Fig. D7; where again, we see that the results agree with PHENOMA.

When we compare the synthetic WDB populations to the sample of verification binaries we note that both BSE and BPASS do have systems in the strain-frequency space of the verification binaries. However, the synthetic populations predict many systems with lower ASD, which are likely to be similar to the verification binaries in mass and period but at greater distances, where they will be more difficult to detect in electromagnetic surveys.

### 3.3 The detectability of the synthetic binary populations

To count the number of WDBs and BHBs that can be detected in our synthetic populations, we count how many are above the LISA sensitivity curve. This is calculated by both PHENOMA and LEGWORK using the method detailed in Appendix C. In this work, we use an integrated observation time of four years with an assumed threshold SNR of 7 to be resolved. We do not impose to have unique systems per frequency bin and discuss the implications of this below.

The SNR of sources versus their frequency are shown in Figs 12 and 13. The number of systems above the SNR of 7 are shown in Table 3, we see that both PHENOMA and LEGWORK predict numbers of detected systems agree closely. In this table, we have also calculated the detected population at four different locations around the model galaxy to determine whether the position of LISA in the galaxy affects the population detected.

Fig. 12 displays the individual SNRs of BHBs in our two populations. BPASS has a greater population of high-frequency BHBs. The result of this is that four sources are predicted to be detected from the BPASS model, compared to only one from the BSE model. We

note the work of Sesana, Lamberts & Petiteau (2020) built on the work of Lamberts et al. (2018), using the same BSE model but made 300 synthetic realisations of the BHB populations and predicted 4.2 detectable BHBs with the same SNR threshold of 7 and observation time. This suggests that in the BSE population, the number of detected systems we employ is 0.5, this might just be because of the distributions of BHBs in this model population rather than an inherent difference between BSE and BPASS. Also, Sesana et al. (2020) found that having one detected BHB was within the 90 per cent confidence interval. A large number of BHB predictions are compiled by Wagg et al. (2022d), with predicted BHB detections from zero to 154, although most studies predict a few BHB detections. Thus, at least in the high-mass regime, BPASS is consistent with other work.

We also plot predictions of individual SNRs of Galactic WDBs in two populations in Fig. 13. The two WDB populations overlap each other most at low frequencies. The BSE model predicts 13 314 WDB sources above the threshold of 7, and BPASS predicts 673 WDBs with an SNR greater than 7. The BPASS WDB population has more systems at lower frequencies. However, BPASS has fewer individually observable WDB systems than that of the BSE WDB population. This significant difference between the BPASS and BSE WDB are discussed below. The BSE results are more similar to other predictions (e.g. Nelemans et al. 2001b) where it is expected that 10 000s of WDBs will be resolved by LISA.

The effect of the choice of the Sun's position on the number of individually resolved systems is also tested, and results are shown in Table 3. The predictions from different choices for the Sun's position in the galaxy do not vary significantly but they allow us to estimate the uncertainty in the detected number of sources due to where LISA is positioned within the Galaxy.

In Table 3, we have also included the number of sources that can be detected if we remove the confusion noise. Give that with BPASS, we predict a factor of 20 fewer WDBs, it is highly likely that the confusion noise from unresolved binaries will be correspondingly lower. Without the confusion noise, the number of observable WDBs increases to 1048. This indicates that for a smaller population of
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The four-year mission SNR for the BHBs in the two model populations versus frequency. We also plot LISA verification binaries for comparison. The SNR threshold marked with dashed line of SNR = 7.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** The four-year mission SNR for WDBs in the two model populations versus frequency. We also plot LISA verification binaries for comparison. Both populations predict systems similar to the verification binaries within the LISA frequency range. The SNR threshold of 7, which is represented by the horizontal dashed line.</figcaption></figure>

WDBs, more individual sources can be detected compared to the populations typically predicted by other studies.

When combined, the results confirm previous findings that WDBs and BHBs will be detected by LISA but the true number is still uncertain. Unfortunately, it appears that significant uncertainties remain in the modelling of the binary evolution itself. Even when choosing different locations within the Galaxy, the number of detected signals does not significantly vary. Although the model of the galaxy itself will have a significant impact on the predictions. Thus, when preparing for LISA and its data analysis, flexibility must be included because until it flies, we cannot be certain about the signals it will detect beyond the known verification binaries.

If we compare to the tabulated results of Wagg et al. (2022c), we find that our predicted BHB population of a few is typical of other predictions (see fig. 12 in Wagg et al. 2022d), although they suggest that up to 154 BHBs (for certain model parameters) could be detected! This is because Wagg et al. (2022d) widely varied the input physics in the COSMIC population synthesis code (an updated user-friendly version of BSE<fnref n="5" />). It is true that this gives an insight into the uncertainty of the binary population synthesis, but some of those parameters would be ruled out by observations. BPASS has been repeatedly validated against numerous different observations (i.e.

<footnote n="5">⁵https://cosmic-popsynth.github.io/COSMIC/</footnote>
<table><caption>**Table 3.** Number of Galactic binary systems with SNR > 7 in the MW-like galaxy predicted by different models, for different locations for the Sun in the galaxy for PHENOMA (top half) and LEGWORK (bottom half).</caption><thead><tr><th rowspan="2">Location<sup>☉</sup></th><th colspan="4">PHENOMA</th></tr><tr><th colspan="2">BHB</th><th colspan="2">WDB</th></tr><tr><th></th><th>BSE</th><th>BPASS</th><th>BSE</th><th>BPASS</th></tr></thead><tbody><tr><td>(8, 0, 0)</td><td>–</td><td>4</td><td>13 296</td><td>682</td></tr><tr><td>(–8, 0, 0)</td><td>1</td><td>5</td><td>13 375</td><td>676</td></tr><tr><td>(0, 8, 0)</td><td>1</td><td>4</td><td>13 187</td><td>665</td></tr><tr><td>(0, –8, 0)</td><td>–</td><td>5</td><td>13 401</td><td>681</td></tr><tr><td>Mean</td><td>0.5±0.6</td><td>4.5±0.6</td><td>13 315±96</td><td>676±8</td></tr><tr><th colspan="5">LEGWORK</th></tr><tr><td>(8, 0, 0)</td><td>–</td><td>4</td><td>13 882</td><td>678</td></tr><tr><td>(–8, 0, 0)</td><td>1</td><td>5</td><td>13 893</td><td>669</td></tr><tr><td>(0, 8, 0)</td><td>1</td><td>4</td><td>13 505</td><td>660</td></tr><tr><td>(0, –8, 0)</td><td>–</td><td>5</td><td>13 959</td><td>674</td></tr><tr><td>Mean</td><td>0.5±0.6</td><td>4.5±0.6</td><td>13 810±206</td><td>670±8</td></tr><tr><td>Combined mean</td><td>0.5±0.5</td><td>4.5±0.5</td><td>13 562±304</td><td>673±8</td></tr><tr><th colspan="5">PHENOMA – no confusion noise</th></tr><tr><th rowspan="2">Location<sup>☉</sup></th><th colspan="2">BHB</th><th colspan="2">WDB</th></tr><tr><th>BSE</th><th>BPASS</th><th>BSE</th><th>BPASS</th></tr><tr><td>(8, 0, 0)</td><td>–</td><td>4</td><td>19 156</td><td>1036</td></tr><tr><td>(–8, 0, 0)</td><td>1</td><td>5</td><td>19 154</td><td>1062</td></tr><tr><td>(0, 8, 0)</td><td>1</td><td>4</td><td>18 867</td><td>1048</td></tr><tr><td>(0, –8, 0)</td><td>–</td><td>5</td><td>19 281</td><td>1045</td></tr><tr><td>Mean</td><td>0.5±0.6</td><td>4.5±0.6</td><td>19 115±175</td><td>1,048±11</td></tr></tbody></table>

Wofford et al. 2016; Eldridge et al. 2017, 2019; Stanway & Eldridge 2018; Xiao et al. 2019; Byrne, Stanway & Eldridge 2021; Massey et al. 2021). While the physics in BPASS can still be refined, the repeated agreement with observed stars gives confidence in results. However, in the low-mass regime of WD evolution, this is one area where further testing must occur, and our results for the WDB population are less certain.

We further remark that when considering the detectability of a system in the LISA frequency band, it is important to consider the SNR and the presence of other binaries within the same frequency bin; especially, the number of systems is high in the lower LISA frequency bins (see Breivik et al. 2020a). In the LISA frequency band, each binary system emits GWs at a specific frequency, assuming the signals are stationary during the four-year LISA mission time), but the LISA instrument has a finite frequency resolution (width of each frequency bin). Even if a binary system has a sufficiently high SNR > 7 to be detected above the noise, it is still crucial to consider the presence of other binaries within the same frequency bin. A full detection pipeline is needed. Suppose another binary system has a similar or identical emission frequency within that bin. In that case, it may be challenging to distinguish the individual signals from each other and from the background noise. In such a scenario, the multiple overlapping signals could interfere with each other, making it difficult to extract the individual signals. This effect is particularly relevant when numerous binaries emit GWs within a specific frequency bin, which is possible given the Galaxy’s large expected population of binary systems. This effect was discussed in more detail in van Zeist et al. (2023). As a result, the number of detectable binaries may be reduced because the overlapping signals hinder the ability to separate them from the noise and the other binaries within the same bin. We need to consider the SNR and the potential interference caused by multiple binaries within the same frequency bin when predicting the

detectability of GWs in the LISA frequency band. This consideration could help provide a more realistic estimation of the number of observable binaries in our research. The size of this can be estimated by comparing our number of detected sources in our study, of 13 500 WDBs to the 12 000 detected by Lamberts et al. (2019) who included this effect, suggesting that approximately 11 per cent fewer systems may be detected.

## 4 DETECTION OF SYNTHETIC BINARY POPULATIONS AND CREATION OF THE MOCK-LISA SIGNAL

Now that we have our population of WDBs and BHBs, our next step is to predict the GW signals that LISA will detect from all the sources at once. The aim is to create the ASD of the Galactic foreground signal from WDBs and BHBs that will be observed over a year. To achieve this, we use the method of Boileau et al. (2021). The full details are described in Appendix A.

We demonstrate in previous sections the detectability by LISA of different binary populations in model galaxies. In this section, we present the results of our modulated signal predictions in time domain. This is to show that we can produce such signals, which will lead to our next data analysis paper. In creating these modulated signals, we include all WDBs and BHBs with a frequency of GWs above $10^{-5}$ Hz and we use the method outlined in Appendix A.

The motion of LISA around the Sun causes a periodic variation in the GW signal detected by LISA, because of the Doppler effect associated with this motion. This variation has a characteristic pattern called a ‘sine-Gaussian’ modulation, from which one can extract information about the binary system, such as sky position, distance, orientation, and orbital parameters, and study the MW galaxy. The potential to extract information about binary systems makes studying the modulation important.

From our simulated binary populations, we produce and predict a modulated GW signal for each binary in the Galaxy according to equation (A12). The time variable changes linearly from 0 s to $3.16 \times 10^7$ s or 1 yr, with a time step of 771 s. Actual LISA data will be recorded with a much smaller sampling frequency. With a 771 s time-step, we require $2^{12}$ data points, which are sufficient to produce the modulation without requiring too much computational time. To produce the linear combination modulation, we sum up all the individual signals in the galaxy to produce the final modulated signals for a particular binary population (using equation A13). We show the mock time domain signals in Figs 14 and 15. Comparing the BHB and WDB populations, the BHB strain signal is about one order of magnitude stronger than that of the WDB population, and this difference is to be expected as the GW strain is proportional to the mass of the binary system.

Both Figs 14 and 15 show modulated signals, as the measured wave amplitude changes over the LISA orbit. Physically, modulation occurs when LISA faces the Galactic centre where the population is concentrated. BPASS and BSE predict that the modulated signal from the WDB is about two orders of magnitude stronger than that of the BHB. Fig. 14 shows the that the maximum strain of the modulated BHB population strain reaches to just below $2 \times 10^{-20}$ for the BPASS model, compared with that of BSE, which is double this value. In comparison, the reverse is true for WDBs. The strain of the modulated WDB population signal for the BPASS model is about twice that predicted by BSE, shown in Fig. 15. For both modulated signals, the amplitude is greatest for the population with the greatest number of sources, rather than being dominated by the highest strain systems.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Modulated GW strain signals of Galactic BHB populations are shown over the duration of a one-year LISA mission. The BPASS line represents the sum of 9298 BPASS Galactic BHBs, and the BSE line represents the sum of 25 335 BSE Galactic BHBs.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Modulated GW strain signals of Galactic WDB populations are shown over the duration of a one-year LISA mission. The BPASS line represents the sum of 58 522 007 BPASS Galactic WDBs, and the BSE line represents the sum of 35 832 029 BSE Galactic WDBs.</figcaption></figure>

The differences between BPASS and BSE are expected as BPASS predicts seven times more Galactic WDB than BSE. While BSE predicts 2.5 times the number of BHBs than BPASS. Because the total strain signal is the sum of all the signals the more numerous the sources, the higher the combined strain. This is also demonstrated in Fig. D8, where we show a total strain including

WDBs with frequencies below $10^{-5}$ Hz which increase the strain further.

These results show that the synthetic signals are sensitively linked to the binary evolution model used. Furthermore, their features are linked in a non-linear way to how many sources are detectable. It is important to note, because analysis pipelines that take this data
and analyse it, must have the flexibility to allow for a wide range of
possible signals.

## 5 DISCUSSION

In this article, we have combined the FIRE high-resolution cosmological simulation and BPASS binary population predictions for GW sources. This has allowed us to predict the number and parameter distributions of unmerged Galactic BHBs and WDBs in a simulated MW-like galaxy. From this, we have simulated the expected modulated LISA signals and identified the expected population of sources that can be resolved with LISA. Our predictions will aid in developing data analysis methods, we show how this will be done in a future paper.

The GW population prediction is sensitive to the choice of binary population synthesis model. We compared our results with the previous work of Lamberts et al. (2018, 2019) who used the same FIRE galaxy simulation but the BSE binary population synthesis code. This allows us for the first time to see how sensitive such predictions are when only the population synthesis code used is varied. We find significantly different populations of WDBs and BHBs, showing that choice of stellar models introduces a significant uncertainty on the predicted compact remnant binary population.

The differences are apparent in all aspects of the results; not just in the total number of sources but in the distribution of masses, frequency, metallicity and the number of systems that can be resolved by LISA. It is key to validate these binary population synthesis codes by comparing to other observations. Multiple observational constraints should be used to counter degeneracies between the input physics and how it affects predictions. We note that the BPASS models have been widely tested against many varied previous observational constraints. However, there were historically relatively few comparisons to lower mass stars relevant for WDBs, though we have begun to explore this regime (e.g. Briel et al. 2022; Byrne et al. 2022).

When we compare the BSE and BPASS, BPASS predicts more detectable BHBs than BSE despite predicting fewer BHBs overall. This is because the mass and frequency distribution of BPASS BHBs extends to higher chirp masses and frequencies. Although we note that Sesana et al. (2020), using more realizations of BSE, did find on average that 4.2 BHBs would be detected. We also note that BPASS reproduces the high-mass end of the GW transient chirp mass distribution (e.g. Briel et al. 2023), which suggests that BPASS predictions for LISA BHBs are closer to what might be observed. However, we note that there are of course other pathways such as dynamical interactions in clusters to form GW transients.

For the WDB population, BPASS predicts more WDBs in total but fewer observable systems than BSE. BPASS has not been validated for such low-mass stars as it has been for high-mass stars. So for WDBs, it is difficult to be certain as to which population is more realistic, at least until more complete catalogues for known WDBs from electromagnetic observations are available. However, we note that BPASS can reproduce the observed type Ia rate and its evolution with redshift (detailed in Eldridge et al. 2019; Briel et al. 2022). However, the majority of Type Ia SNe would then be via the single-degenerate rather than double-degenerate channel.

Binary evolution is uncertain and the different populations found for the WDBs here show that different input physics can lead to order of magnitude differences in the predicted number of LISA sources. Further work to understand these differences will require the comparison of a large number of the extant binary population

synthesis codes as well as varied observations to constrain them until LISA launches and takes its first observations.

This study is one of the first to directly compare the difference in outputs between a binary stellar evolution models based on BSE rapid stellar evolution (Hurley et al. 2002) and that from the BPASS detailed stellar evolution models (Stanway & Eldridge 2018). It is an interesting comparison because the equations that represent the stellar evolution for BSE (Hurley et al. 2000) were based on detailed stellar evolution models calculated by the Cambridge STARS code originally described in Eggleton (1971) with the last paper in common between the two being Pols et al. (1995).

Given this fact, one might expect that the results from the two codes would be identical for stellar populations. But, they are clearly not. The reason for the differences is complex. While for single stars, the two codes must create similar results due to the same stellar physics being used. There are two differences between the BPASS version of the Cambridge STARS code and that of Pols et al. (1995), BPASS uses newer opacity tables (Eldridge & Tout 2004) and an updated stellar wind mass-loss prescription (Eldridge et al. 2017). The differences from this however only impact on the latest stages of stellar evolution and the BHBs from massive stars. Thus, differences are thus likely to arise from the different ways binary interactions are dealt with.

The binary physics model of Hurley et al. (2002) was used as a guide for the physics of BPASS but several modifications had to be made because of the nature of how detailed stellar evolution models work. For example, it is not possible to just remove the hydrogen envelope in a detailed model, it must be done over a number of time-steps. Thus, the CEE prescription is inherently different. The implementation of the CEE within BPASS also assumes angular momentum conservation rather than the typical conservation of energy that is employed in BSE (see the supplementary data of Stevance et al. 2023, for further discussion on the different CEE prescription). Thus, one possible explanation for fewer tight WDBs in BPASS is that the CEE model is not producing enough tight binaries. We have made some recent exploratory models that indicate that tighter WDBs can be made if a BSE-like CEE prescription is implemented within BPASS, so this is the most likely cause of the difference between BPASS and BSE WDB populations. However, creating a full population to test this would require to recalculate a new series of BPASS stellar evolution models with this different CEE scheme. This would be highly computationally expensive but is clearly something we must further investigate.

Another important source of the differences between the BSE and BPASS results is the relative rates of SMTs. In BSE, mass transfer stability is determined by analytic expressions comparing the evolutionary state of the stars and their mass ratio (Hurley et al. 2002). In BPASS, when the donor star fills its Roche lobe, mass is lost from the donor and the donor star then either reaches an equilibrium so that the stellar radius is stable, or the stellar radii continue to grow so that CEE occurs. Thus, in BPASS, how the star responds to mass-loss determines what happens to the binary. We find, especially at lower masses, MT is much more stable than would be assumed in BSE, leading to a much lower occurrence of CEE. Thus, even though the CEE models are different, it is the fact that SMT is more prevalent that could be the reason for much of the difference. The clearest indication of this difference is the case shown in Fig. 5. Here, the fact so many of the secondary WDs are more massive than the primary WDs is due to efficient MT during the first occurrence of MT within the binary. Further implications of the prevalence of SMT in BPASS are discussed in Briel et al. (2023).
A final difference between BSE and BPASS comes from a more subtle result of the use of detailed stellar evolution models. MT can occur at almost any phase of a star's evolution. How the core mass and core convective zones react to changes in a star's mass can be accurately followed in a detailed stellar evolution code. MT also ends when a star shrinks within its Roche lobe rather than assuming a star must lose its whole hydrogen envelope. These differences in how the stellar models evolve could be the reason behind the different structures in the chirp mass–GW frequency distributions of Fig. 7. Especially, at the lower chirp masses, where BSE does not allow some WDBs to exist while they are predicted in the BPASS results.

All of these factors lead to the different distribution of WDBs from BPASS that we find compared to that from Lamberts et al. (2018, 2019). In particular, we find that there are fewer CEE events and mergers from stars with initial masses below approximately 5 M☉. These low-mass stars have much slower evolutionary time-scales and with our method of determining MT stability this leads to mass transfer being much more stable as the stars do not continue to grow in radius when they fill their Roche lobes. This increased MT stability means they avoid CEE, leading to fewer very close orbits and fewer binaries merging.

Further analysis of the source of the difference between the BPASS and BSE binary evolution algorithms is beyond the scope of this work. However, we can draw a firm conclusion that the treatment of physical processes within binary population synthesis codes does have a significant impact on the predicted WDB and BBH populations. This indicates that the future observations by LISA will play a important role in constraining our understanding of those physical processes.

## 6 CONCLUSION

We summarize our most important findings in this paper as follows:

(i) Our simulations predict a population of BHBs in the MW. The BPASS population is smaller than that of Lamberts et al. (2018). However, the number of detectable BHBs, with an SNR greater than 7, is greater. We expect to detect $4.5 \pm 0.5$ from BPASS predictions compared to $0.5 \pm 0.5$ from BSE. Our prediction is similar to Sesana et al. (2020), whom predicted 4.2 detectable BHB systems from using a large number of BSE model WDB populations. Our values also agree with the simulations collated by Wagg et al. (2022d). With the range of values found by other works ranging from 0 to 154.

(ii) Our results confirm that there will be a significant overall population of WDBs in the MW galaxy, greater than the previous study (Lamberts et al. 2019) found. The predicted distribution of masses and separations of these WDBs is consistent with observations of LISA verification binaries and the extrapolated double-degenerate WDB population. However, due to the verification binaries being an inhomogenous sample, it is difficult to compare our predicted number of such binaries beyond 10 kpc.

(iii) We predict $673 \pm 8$ LISA resolvable Galactic WDBs from BPASS compared with $13562 \pm 304$ from BSE that have an SNR greater than 7 and would be able to be resolved by LISA (see Table 3).

(iv) With the smaller number of WDBs from BPASS, it is likely that the confusion noise from undetected WDBs will be substantially weaker. If we assume this means that we can ignore the confusion noise BPASS then predicts up to $1048 \pm 11$ WDBs could be detected.

(v) We find our prediction for WDB systems with SNR greater than 7 is only weakly dependent on the tool used to determine detectability. With predicted numbers varying a few per cent between

PHENOMA and LEGWORK. We see no difference in our BHB numbers but there are very few objects detected in this population.

(vi) We find the choice of binary population synthesis codes introduces significant differences in the predicted compact remnant binary population. This is starkly revealed in Figs 10 and 11. This is likely due to differences in how the binary interactions are dealt with in BSE and BPASS, as well as due to the use of rapid stellar evolution modelling in BSE compared to the detailed modelling of the stellar evolution in BPASS. Further comparison to low-mass binary stars with WDs will be key to revealing which synthetic binary population is closer to what LISA will be expected to observe. We also need to evaluate the different pathways that lead to WDBs and see what impact extra physics such as magnetic-wind braking by white dwarfs may have on our predicted populations (e.g. Van, Ivanova & Heinke 2019). We do find that the total number of predicted WDBs by BPASS is similar to what is expected based on the density of WDBs in the Solar neighbourhood.

Finally, we note that the simulation-based method we have used provides a powerful tool for predicting the population of BHBs and WDBs in an MW mass-like galaxy. Our predicted populations are consistent with other observational constraints, but the low-mass models need more validation. This is difficult because binary population synthesis codes have now been extant for decades, and providing strong constraints of the evolution of binary star populations is still somewhat uncertain. In this work, we have demonstrated that LISA observation will provide an important constraint on binary populations. As detailed in the LISA Mission L3 proposal (see details in recent report by Colpi et al. 2024), the data processing for LISA involve several distinct levels of analysis to extract meaningful information from the raw data. The insights and mock data we have created here may have use as future LISA data flows are designed (at levels 2 and 3).

## ACKNOWLEDGEMENTS

We thank Max Briel for useful discussions on MT in BPASS models. We gratefully acknowledge support by the Marsden Fund Council from New Zealand Government funding, managed by Royal Society Te Apārangi, the University of Auckland and the University of Warwick for their continuing support. BPASS is enabled by the resources of the NeSI Pan Cluster. New Zealand's national facilities are provided by the NZ eScience Infrastructure and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation & Employment's Research Infrastructure programme. URL: https://www.nesi.org.nz. AL acknowledges support by the ANR COSMERGE project, grant ANR-20-CE31-001 and the 'Programme National des Hautes Energies' (PNHE) of CNRS/INSU co-funded by CEA and CNES. GB acknowledges support from the Centre National D'Études Spatiales (CNES).

## DATA AVAILABILITY

The BPASS Galactic binary populations used in this article will be shared on reasonable request to the corresponding author. For more a general BPASS output regarding binary evolution, please visit https://bpass.auckland.ac.nz or http://warwick.ac.uk/bpass. The PHENOMA model used in this study can be downloaded from the GitHub reference in Robson et al. (2019).

The LEGWORK model used in this study can be downloaded from GitHub, and the detail of the software is outlined in Wagg et al. (2022a).
## REFERENCES

Abbott R. et al., 2021, Phys. Rev. X, 13, 041039

Ajith P. et al., 2007, Class. Quantum Gravity, 24, S689

Amaro-Seoane P. et al., 2023, Living Rev. Relativ., 26, 2

Andrews J. J., Breivik K., Pankow C., D'Orazio D. J., Safarzadeh M., 2020, ApJ, 892, L9

Belczynski K., Bulik T., Fryer C. L., 2012, preprint (arXiv:1208.2422)

Belczynski K., Bulik T., Mandel I., Sathyaprakash B. S., Zdziarski A. A., Mikołajewska J., 2013, ApJ, 764, 96

Boileau G., Lamberts A., Christensen N., Cornish N. J., Meyer R., 2021, MNRAS, 508, 803

Breivik K. et al., 2020a, ApJ, 898, 71

Breivik K., Mingarelli C. M. F., Larson S. L., 2020b, ApJ, 901, 4

Briel M. M., Eldridge J. J., Stanway E. R., Stevance H. F., Chrimes A. A., 2022, MNRAS, 514, 1315

Briel M. M., Stevance H. F., Eldridge J. J., 2023, MNRAS, 520, 5724

Burdge K. B. et al., 2020, ApJ, 905, 32

Byrne C. M., Stanway E. R., Eldridge J. J., 2021, MNRAS, 507, 621

Byrne C. M., Stanway E. R., Eldridge J. J., McSwiney L., Townsend O. T., 2022, MNRAS, 512, 5329

Cantiello M., Yoon S. C., Langer N., Livio M., 2007, A&A, 465, L29

Casares J., Jonker P. G., Israelian G., 2017, in Alsabti A. W., Murdin P., eds, Handbook of Supernovae. Springer, Berlin, p. 1499

Colpi M. et al., 2024, preprint (arXiv:2402.07571)

Cornish N. J., 2001, Phys. Rev. D, 65, 022004

Cornish N. J., Larson S. L., 2001, Class. Quantum Gravity, 18, 3473

Cornish N. J., Littenberg T. B., 2007, Phys. Rev. D, 76, 083006

Cornish N. J., Littenberg T. B., Bécsy B., Chatzioioannou K., Clark J. A., Ghonge S., Millhouse M., 2021, Phys. Rev. D, 103, 044006

Eggleton P. P., 1971, MNRAS, 151, 351

Eldridge J. J., Stanway E. R., 2022, ARA&A, 60, 455

Eldridge J. J., Tout C. A., 2004, MNRAS, 348, 201

Eldridge J. J., Izzard R. G., Tout C. A., 2008, MNRAS, 384, 1109

Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, PASA, 34, e058

Eldridge J. J., Stanway E. R., Tang P. N., 2019, MNRAS, 482, 870

Farmer A. J., Phinney E. S., 2003, MNRAS, 346, 1197

Finch E. et al., 2023, MNRAS, 522, 5358

Gaia Collaboration, 2018, A&A, 616, A11

Gallegos-Garcia M., Berry C. P. L., Marchant P., Kalogera V., 2021, ApJ, 922, 110

Ghodla S., Eldridge J. J., 2023, MNRAS, 523, 1711

Ghodla S., van Zeist W. G. J., Eldridge J. J., Stevance H. F., Stanway E. R., 2022, MNRAS, 511, 1201

Hannam M., 2014, Gen. Relativ. Gravit., 46, 1767

Hils D., Bender P. L., Webbink R. F., 1990, ApJ, 360, 75

Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974

Holberg J. B., Oswalt T. D., Sion E. M., McCook G. P., 2016, MNRAS, 462, 2295

Hopkins P. F., 2015, MNRAS, 450, 53

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543

Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897

Karnesis N., Babak S., Pieroni M., Cornish N., Littenberg T., 2021, Phys. Rev. D, 104, 043019

Klencki J., Nelemans G., Istrate A. G., Chruslinska M., 2021, A&A, 645, A54

Kobayashi C. et al., 2023, ApJ, 943, L12

Korol V., Rossi E. M., Groot P. J., Nelemans G., Toonen S., Brown A. G. A., 2017, MNRAS, 470, 1894

Korol V., Hallakoun N., Toonen S., Karnesis N., 2022, MNRAS, 511, 5936

Kremer K., Breivik K., Larson S. L., Kalogera V., 2017, ApJ, 846, 95

Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545

Kupfer T. et al., 2018, MNRAS, 480, 302

Kupfer T. et al., 2023, ApJ, 963, 14

Lamberts A. et al., 2018, MNRAS, 480, 2704

Lamberts A., Blunt S., Littenberg T. B., Garrison-Kimmel S., Kupfer T., Sanderson R. E., 2019, MNRAS, 490, 5888

Laplace E., Justham S., Renzo M., Göteborg Y., Farmer R., Vartanyan D., de Mink S. E., 2021, A&A, 656, A58

Lau M. Y. M., Mandel I., Vigna-Gómez A., Neijssel C. J., Stevenson S., Sesana A., 2020, MNRAS, 492, 3061

Long A. S. et al., 2020, ApJ, 898, 133

Marchant P., Pappas K. M. W., Gallegos-Garcia M., Berry C. P. L., Taam R. E., Kalogera V., Podsiałowski P., 2021, A&A, 650, A107

Massey P., Neugent K. F., Dorn-Wallenstein T. Z., Eldridge J. J., Stanway E. R., Levesque E. M., 2021, ApJ, 922, 177

Moe M., Di Stefano R., 2017, ApJS, 230, 15

Nelemans G., Yungelson L. R., Portegies Zwart S. F., Verbunt F., 2001a, A&A, 365, 491

Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001b, A&A, 375, 890

Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2004, MNRAS, 349, 181

Nissanke S., Vallisneri M., Nelemans G., Prince T. A., 2012, ApJ, 758, 131

Patton R. A., Sukhbold T., Eldridge J. J., 2022, MNRAS, 511, 903

Pelisoli I. et al., 2021, Nat. Astron., 5, 1052

Peters P. C., 1964, Phys. Rev., 136, 1224

Pols O. R., Tout C. A., Eggleton P. P., Han Z., 1995, MNRAS, 274, 964

Prince T. A., Tinto M., Larson S. L., Armstrong J. W., 2002, Phys. Rev. D, 66, 122002

Rappaport S., Verbunt F., Joss P. C., 1983, ApJ, 275, 713

Riley J. et al., 2022, J. Open Source Softw., 7, 3838

Robson T., Cornish N. J., Liu C., 2019, Class. Quantum Gravity, 36, 105011

Rosswog S., Kasen D., Guillochon J., Ramirez-Ruiz E., 2009, ApJ, 705, L128

Ruiter A. J., Belczynski K., Benacquista M., Larson S. L., Williams G., 2010, ApJ, 717, 1006

Schmidt P., Ohme F., Hannam M., 2015, Phys. Rev. D, 91, 024043

Sesana A., Lamberts A., Petiteau A., 2020, MNRAS, 494, L75

Seto N., 2002, MNRAS, 333, 469

Seto N., 2022, Phys. Rev. Lett., 128, 041101

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stevance H. F., Eldridge J. J., Stanway E. R., Lyman J., McLeod A. F., Levan A. J., 2023, Nat. Astron., 7, 444

Stroeer A., Vecchio A., 2006, Class. Quantum Gravity, 23, S809

Thiele S., Breivik K., Sanderson R. E., Luger R., 2023, ApJ, 945, 162

Van K. X., Ivanova N., Heinke C. O., 2019, MNRAS, 483, 5595

van Son L. A. C. et al., 2020, ApJ, 897, 100

van Zeist W. G. J., Eldridge J. J., Tang P. N., 2023, MNRAS, 524, 2836

van den Heuvel E. P. J., 2019, IAU Symp., 346, 1

Wagg T., Breivik K., de Mink S., 2022a, J. Open Source Softw., 7, 3998

Wagg T., Breivik K., de Mink S. E., 2022b, ApJS, 260, 52

Wagg T., Breivik K., de Mink S. E., 2022c, ApJS, 260, 52

Wagg T., Broekgaarden F. S., de Mink S. E., Frankel N., van Son L. A. C., Justham S., 2022d, ApJ, 937, 118

Wetzel A., Garrison-Kimmel S., 2020, Gizmo Analysis: Read and analyze Gizmo simulations, Astrophysics Source Code Library. preprint (ascl:2002.015)

Wetzel A. et al., 2023, ApJS, 265, 44

Wofford A. et al., 2016, MNRAS, 457, 4296

Xiao L., Galbany L., Eldridge J. J., Stanway E. R., 2019, MNRAS, 482, 384

Yoon S. C., Langer N., 2005, A&A, 443, 643

## APPENDIX A: GW WAVEFORM CALCULATION

For the benefit of readers interested in exploring the data and results of this study in greater depth, this appendix presents mathematical formulations, alternative versions of several key plots and supplementary plots included in the main body of the paper. These alternate plots offer additional insights and perspectives that enhance the reader's understanding of the research while providing a transparent and comprehensive view of the study's findings.
The first step in calculating the GW signal from a compact binary is to predict the waveform from these objects. In an MW-like galaxy, the distance between a binary system of a Cartesian coordinate (X, Y, Z) and LISA ($X_\odot, Y_\odot, Z_\odot$) = (8, 0, 0) in kpc is

$$R = \sqrt{(X - X_\odot)^2 + (Y - Y_\odot)^2 + (Z - Z_\odot)^2}. \qquad (A1)$$

The strain of the GW signal is proportional to the amplitude of the signal; and the GW amplitude depends on the mass of the binary and the distance between the source and the detector. Similar to electromagnetic radiation, GWs also come with two polarizations, namely plus (orthogonal) and cross (diagonal). The amplitude of the plus and cross-polarized GW emitted by a binary is (Boileau et al. 2021)

$$A_+(M_1, M_2, R, f, t) = \frac{2G^2 M_1 M_2}{c^4 R} \left( \frac{(\pi f)^2}{G(M_1+M_2)} \right)^{1/3} \times (1 + \cos^2(t)) \qquad (A2)$$

$$A_\times(M_1, M_2, R, f, t) = -\frac{2G^2 M_1 M_2}{c^4 R} \left( \frac{(\pi f)^2}{G(M_1+M_2)} \right)^{1/3} 2 \cos(t). \qquad (A3)$$

where $M_1$ and $M_2$ are primary and secondary masses, in solar mass M$_\odot$, of the binary. G is the Gravitational constant. f is the GW frequency, which is twice the binary's orbital frequency for circular orbits. The inclination t is the angle between the orbital plane of the binary and an observer. A 0° inclination means the GW is travelling face-on towards LISA and has maximum strain signal. Assuming maximum inclination then a factor of 4/5 needs to be included in the equation (A3) if one wants to consider the average strain over all inclination angles (Robson et al. 2019). Studies such as Finch et al. (2023) have shown from their models that different inclinations can affect amplitudes estimation. Another method is to draw $\cos(t)$ from a uniform distribution as there is no way to measure individual binary star inclination (Boileau et al. 2021). In our models, we assume $\cos(t)$ is randomly drawn from a uniform distribution between -1 and 1.

The time-dependent GW strain $h(t)$ is given by the linear combination of two polarizations of the waveform,

$$h(t) = h(t)_+ e_+ + h(t)_\times e_\times, \qquad (A4)$$

where $e_+$ and $e_\times$ are the two polarization tensors, which are given by

$$\mathbf{e}_{+} = \mathbf{E} \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{pmatrix} \mathbf{E}^{T} \quad \mathbf{e}_{\times} = E \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \mathbf{E}^{T}. \qquad (A5)$$

The polarization coordinate matrix E is given as

$$\mathbf{E} = \begin{pmatrix} \lambda_s \psi_c - \lambda_c \beta_s \psi_s & -\lambda_s \psi_s - \lambda_c \beta_s \psi_c & -\lambda_c \beta_s \\ -\lambda_c \psi_c - \lambda_s \beta_s \psi_s & \lambda_c \psi_s - \lambda_s \beta_s \psi_c & -\lambda_s \beta_c \\ \beta_c \psi_s & \beta_c \psi_c & -\beta_s \end{pmatrix}, \qquad (A6)$$

where the subscripts s and c are sine and cosine operations, respectively, the Galactic coordinates are latitude $\beta$ and longitude $\lambda$, and $\psi$ is the rotation around the direction of GW propagation; as the effect of varying the $\psi$ values is not the focus of this paper, to simplify computation, we have set $\psi$ to 0.

Many Galactic binaries have a frequency within the LISA frequency band and do not merge; however, some tight orbit systems show frequency variation. The computation for systems with frequency variation is more complex. For simplicity, ignoring systems that have frequency variation within the four-year LISA mission

<figure><img src="image_17.png" /><figcaption>**Figure A1.** LISA coordinate system ($\hat{u}, \hat{v}$) used to evaluate the transfer function which is discussed in Cornish & Larson (2001). LISA spacecrafts are located at the triangle's vertices, and all three arms are equal length.</figcaption></figure>

time, here we assume a monochromatic signal. We then have

$$\begin{pmatrix} h_+(t) \\ h_\times(t) \end{pmatrix} = \begin{pmatrix} A_+(f) \cos(2\pi f t + \phi_0) \\ A_\times(f) \sin(2\pi f t + \phi_0) \end{pmatrix}. \qquad (A7)$$

The change in frequency or drift frequency $\dot{f}$ is small over the LISA mission time and ignored in our computational model; though studies, such as Seto (2022), suggest that it is possible to measure the drift frequency of Galactic ultracompact binaries. However, for most systems, $\dot{f}$ is small unless f is large. Lamberts et al. (2019) suggest that in a four-year mission that 750 sources may have detectable $\dot{f}$ and thus, have measurable chirp masses. Importantly in our analysis, it will not change the number of detected systems because for systems when $\dot{f}$ is high, the sources are well above the sensitivity curve. This is consistent with other studies such as, Breivik, Mingarelli & Larson (2020b) and Korol et al. (2022). To simulate the signal, for simplicity, we consider the phase $\phi_0$ is drawn from a uniform distribution between 0 and $2\pi$.

Depending on the location of the GW source in the sky, the signals detected by LISA are modulated. To calculate the response of the LISA detector arms, we use (Cornish & Littenberg 2007)

$$\begin{aligned} \mathcal{H}_+(t) &= A_+(f) \cos(2\pi f t + \phi_0) e_+ : D \\ \mathcal{H}_\times(t) &= A_\times(f) \cos(2\pi f t + \phi_0) e_\times : D, \end{aligned} \qquad (A8)$$

where the geometry dependent D is the one-arm detector tensor

$$D = \frac{1}{2} (\hat{u} \otimes \hat{u} - \hat{v} \otimes \hat{v}) \qquad (A9)$$

and $\hat{u}$ and $\hat{v}$ are coordinate systems of the unit vectors along each arm of LISA spacecraft, defined as

$$\hat{u} = (1/20\sqrt{3}/2) \quad \hat{v} = (-1/2\sqrt{3}/20). \qquad (A10)$$

The sky-location-dependent polarized detector response functions outlined in Cornish (2001) are

$$\begin{aligned} F_+ &= -\frac{\sqrt{3}}{4} (1 + \cos(\theta))^2 \sin(2\phi) \\ F_\times &= -\frac{\sqrt{3}}{2} \cos(\theta) \cos(2\phi), \end{aligned} \qquad (A11)$$

where the source location $\phi$ and $\theta$ are related to the angles of the LISA transfer function as discussed in Cornish & Larson (2001). They are related via $\hat{u} \cdot \hat{z} = \sin(\phi + \pi/6) \sin\theta$ and $\hat{v} \cdot \hat{z} = \sin(\phi - \pi/6) \sin\theta$, where $\hat{z}$ is the z direction in the Cartesian coordinate
system. The LISA coordinate system is shown in Fig. A1. We assume that $\theta$ and $\phi$ change linearly and continuously with time in our computation.

Finally, we multiply the polarized strain and response function to generate a modulated signal from one source. To generate a Galactic binary signal from the whole population with N systems, we sum up all linear combinations of the two polarizations of the sources to get

$$s(t) = \sum_{i=1}^{N} \sum_{A=+,-} h_{A,i} (f_{\text{orb},i}, M1_i, M2_i, X_i, Y_i, Z_i, t) \times F_A(\theta, \phi, t) D(\theta, \phi, t) \quad \text{(A12)}$$

or

$$s(t) = \sum_{i=1}^{N} \sum_{A=+,-} \mathcal{H}_{A,i}(t) \times F_A(\theta, \phi, t). \quad \text{(A13)}$$

The full waveform will be used to calculate the modulated signal in later sections, while the predicted strain can be compared to LISA's sensitivity to determine its detectability.

## APPENDIX B: LISA SENSITIVITY CURVE

One way to show LISA's ability to detect GW sources is to compare the amplitude spectral density of the GW signals with the LISA full sensitivity curve $S(f)$ in the frequency domain, which is given in Robson et al. (2019) as

$$S(f) = S_n(f) + S_c(f). \quad \text{(B1)}$$

A full sensitivity curve is the combined PSD of the LISA instrumental noise $S_n$ and the Galactic foreground confusion noise $S_c$. $S_n$ can be defined as

$$S_n(f) = \frac{P_n(f)}{R(f)}, \quad \text{(B2)}$$

where $P_n(f)$ is the power spectral density (PSD) of the detector noise and $R(f)$ is the sky-averaged detector response function in the frequency domain. $R(f)$ is defined as

$$R(f) = \frac{3}{10} \frac{1}{(1 + 0.6(f/f_*)^2)}, \quad \text{(B3)}$$

where $f_* = c/(2\pi L)$ is the transfer frequency (Prince et al. 2002) and L is the LISA arm length of 2.5 million km. It is assumed that all three arms are equal in length. $P_n(f)$ is formed from the approximated single-link optical metrology noise $P_{\text{oms}}$,

$$P_{\text{oms}} = (1.5 \times 10^{-11} \text{ m})^2 \text{ Hz}^{-1}, \quad \text{(B4)}$$

and the single-test-mass acceleration noise $P_{\text{acc}}$,

$$P_{\text{acc}} = (3 \times 10^{-15} \text{ m s}^{-2})^2 \left(1 + \left(\frac{0.4 \text{mHz}}{f}\right)^2\right) \text{Hz}^{-1}. \quad \text{(B5)}$$

The above two types of instrument noises are then combined to form the total LISA instrument noise in a Michelson-style LISA data channel as

$$P_n(f) = \frac{P_{\text{oms}}}{L^2} + 2 \left(1 + \cos^2(f/f_*)\right) \frac{P_{\text{acc}}}{(2\pi f)^2 L^2}. \quad \text{(B6)}$$

To formulate the full LISA sensitivity curve, we combine $S_n$ and the Galactic confusion noise $S_c$ from Galactic WDBs. The signal is modulated by the change in LISA's antenna response as the detector rotates over the LISA mission time. An averaged sensitivity curve approximation is made for our calculation; for an alternative

approximation see Karnesis et al. (2021). For plotting purposes, we define the confusion noise as Robson et al. (2019),

$$S_c(f) = A f^{-7/3} e^{-f^{\alpha} + \beta f \sin(\kappa f)} [1 + \tanh(\gamma (f_k - f))] \text{Hz}^{-1}, \quad \text{(B7)}$$

where the amplitude A is $9 \times 10^{-45}$ and the parameters of analytic fit for a four-year LISA mission time are $\alpha = 0.138$, $\beta = -221$, $\kappa = 521$, $\gamma = 1680$ and $f_k = 0.00113$ (Robson et al. 2019). It is important to appreciate that in our work, as we are modelling the WDB signal, the confusion noise from these binaries is model dependent, as we show later. However, we can show that individual binaries will be above this model limit. Each GW frequency (twice the orbital frequency) in our models is unique to a particular binary, and we assume it does not change over the LISA mission time (very ideal cases). A strain signal from a Galactic binary, including e.g. any of the verification binaries, is calculated as the one-sided, angle-averaged PSD $S_h(f)$

$$S_h(f) = \frac{A^2(f)}{2T}, \quad \text{(B8)}$$

where $A(f)$ is the amplitude of the wave and T is the LISA observation time frequency domain, in our case T is four years. $A(f)$ is related to distance between LISA and the signal source, the frequency and the $M_c$ of the binary system. In later section, we defined the polarized strain which we will use to approximate signals in the time domain. For plotting in the frequency domain, we defined and computed the $A(f)$ using equations (20–23) in Robson et al. (2019).

The amplitude spectral density (ASD) is widely used by the ground-based detector community, for LISA we define the polarization-averaged ASD as

$$ASD(f) = \sqrt{S_h(f)}. \quad \text{(B9)}$$

Each individual ASD of a Galactic binary source is compared with the LISA sensitivity and known verification binaries in this section. To compute and plot $S_n$ and $S_h$ for the Galactic binaries and verification binaries, we use the equations (1–9) within the waveform model PHENOMA (Ajith et al. 2007) and additional PYTHON modules (Robson et al. 2019). There are more generic and accurate models, such as the latest PHENOMB model (Hannam 2014; Schmidt et al. 2015), which calculates spin-precession effects. PHENOMA, however, is sufficient for calculating strains and estimating signal-to-noise ratios (SNRs) in our study.

## APPENDIX C: CALCULATION OF SNR

The SNR detected by the LISA mission is affected by several factors. Some of the critical factors are as follows: First, the SNR is directly proportional to the amplitude of the GW signal (for systems with known frequencies). Therefore, stronger signals produce higher SNR. Visually, the SNR of each binary is given by the ratio of the height of the sensitivity curve and the height of the corresponding ASD in Figs 10 and 11. The LISA mission is sensitive to GWs in the frequency range of 0.1 to 1 mHz. The SNR is higher for signals that fall within this frequency range. Also, the longer the signal persists in the detector's bandwidth, the higher the SNR will be; hence, a longer mission time increases the value of SNR. The SNR is also inversely proportional to the noise level of the detector. The lower the noise level, the higher the SNR will be. Finally, the distance and sky-location of the source also affect the SNR. The closer the source, the higher the SNR will be. Within our models, we can estimate the expected SNR for a given GW signal in the LISA mission.
A binary is detected if the SNR integrated over its observation time (in our case, we only consider the observation time of four years) is larger than an assumed threshold of 7. As the sky-location can affect the SNR for a source (Seto 2022), without including the effect of sky-location, we calculate the SNR of a quasi-circular non-spinning comparable-mass binaries with only a dominant quadrupole harmonic (Robson et al. 2019) of every binary as

$$ \overline{\rho^2} = \frac{16}{5} \int_{f_{\text{in}}}^{f_{\text{in}}+\Delta f} \frac{A^2(f)}{S_n(f)} df \approx \frac{16 \Delta f A^2(f_{\text{in}})}{5 S_n(f_{\text{in}})}, \quad (C1) $$

where $f_{in}$ is the initial GW frequency of the binary. Since we assume sources are monochromatic signals over the LISA mission time, the change in frequency $\Delta f$ can be ignored. We can simplify the SNR to form an approximation as

$$ \overline{\rho^2} \approx \frac{h_{\text{GB}}^2(f_{\text{in}})}{S_n(f_{\text{in}})}, \quad (C2) $$

where $h_{GB}$ is the strain spectral density in Hz$^{-1/2}$, which can be expressed as

$$ h_{\text{GB}} = \frac{8T^{1/2} (G M_c / c^3)^{5/3} \pi^{2/3} f_{\text{in}}^{2/3}}{5^{1/2} (R/c)}, \quad (C3) $$

where R is the distance between the GW source and the Sun.

<figure><img src="image_18.png" /><figcaption>Figure D1. The small-scale 2D number density map of BHB populations in the simulated galaxy from the two model predictions viewed face on. BHBs are present in Galactic bulge and discs in both model predictions. The population density is the highest at the centre of both galaxies.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure D2. The small-scale 2D number density map of WDB populations in the simulated galaxy from the two model predictions viewed face on. WDBs are present in the Galactic bulge and discs in both model predictions. The population density is the highest at the centre of both galaxies.</figcaption></figure>

For more complicated sources with known sky-location and inclination, Robson et al. (2019) suggest a new expression in their equation (36), which calculates SNRs as for the sky-averaged case. Constant SNR scales with $M_c$ to $M_c^{5/3}$ in a given LISA's volume sensitivity, and an SNR threshold of 7 which is computed assuming a four-year LISA mission duration with 100 per cent duty cycle and for (quasi-)stationary equal-mass circular binaries (Amaro-Seoane et al. 2023). LISA will see Sources with SNRs > 7. Extreme mass ratio (EMR) binaries, which may be highly eccentric, require a more involved formalization. Even for binaries with spins, the SNR we compute should be accurate to within a factor of two or so for spinning systems according to Robson et al. (2019).

## APPENDIX D: EXTRA PROPERTIES OF THE BPASS BINARY POPULATIONS

Figs D1 and D2 show small-scale 2D density maps of the comparison between BPASS and BSE for Galactic BHB and WDB populations. The BPASS BHB population is less spread out in the MW than the BSE BHB population, and on the other hand, the BPASS WDB populations are similarly distributed in the MW galaxy.

Most of the BPASS BHBs with $M_c$ greater than 40 M$_{\odot}$ are outside the frequency range of LISA, shown in Fig. D3, BPASS predicts a total number of 495,586 BHBs in the MW galaxy and
<figure><img src="image_20.png" /><figcaption>**Figure D3.** The distribution of BHB by $M_c M_\odot$. The left panel includes the distribution of BHBs over the full frequency range of $10^{-10}-0.1$ Hz, the right panel contain the frequency range ($10^{-5}$ Hz–0.1 Hz).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure D4.** The number distribution of WDB by $M_c M_\odot$. In both panels, the BSE results only include the frequency range $10^{-5}-0.1$ Hz. However, on the left panel, the BPASS results include the full frequency range of $10^{-10}-0.1$ Hz, the right panel contain the frequency range $10^{-5}-0.1$ Hz.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure D5.** The BPASS BHBs and WDBs distribution in $M_c$ versus age distribution.</figcaption></figure>

only 9,298 (1.9 per cent) of the total BHB population are within the LISA frequency band, compared with Lamberts et al. (2018) who found 4.2 per cent within the band. On the other hand, the BPASS predicts a total number of 272 789 350 WDBs and 58 522 007 of these are within the LISA frequency band, nearly 21.5 per cent.

Shown in Fig. D4 range, the population of these ‘non-LISA’ WDBs distributes evenly across the $M_c$, and they contribute to the confusion noise.

We take another look at the BPASS BHB and WDB populations in the LISA frequency band in Fig. D5; where most of the BHBs
<figure><img src="image_23.png" /><figcaption>**Figure D6.** Same as in Fig. 10 but now using LEGWORK. The extra sources above the LISA sensitivity curve may be caused by different inclination assumptions. The inclination is not taken into account in LEGWORK. However, in PHENOMA, we have generated random inclinations to calculate the waveform.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure D7.** Same as in Fig. 11 but now using LEGWORK. The extra sources above the LISA sensitivity curve may be caused by different inclination assumptions. The inclination is not taken into account in LEGWORK. However, in PHENOMA, we have generated random inclinations to calculate the waveform.</figcaption></figure>

have $M_c$ less than 10 M$_{\odot}$ and are about 10 Gyr old. Similarly, older WDBs with $M_c$ below 0.6 M$_{\odot}$ are more numerous in the WDB population.

We provide an alternative frequency versus ASD plot using LEGWORK in Fig. D6. Compared with Figs 10 and D7, we notice a shift in the scale of the ASD. This is because of the difference in the ASD formulation. The alternative ASD defined in LEGWORK is

$$ASD = h_2 \sqrt{T} \qquad (D1)$$

where $h_2$ is the strain amplitude of the source for the second orbital frequency harmonic (Wagg et al. 2022c; Thiele et al. 2023). Both BSE and BPASS BHB populations are similar with a small

fraction of the population above the LISA sensitivity curve at lower frequencies. Unlike the results from the PHENOMA, here the high-frequency WDBs in both BSE and BPASS models have ASDs of about $10^{-17}Hz^{-1/2}$.

We present an alternative WDB population f versus strain plot using LEGWORK in Fig. D7. Similar to the PHEOMA results, the BSE WDB population have similar strains and more systems at a higher frequency than BPASS.

Modulated signal observed by LISA is the superposition of many binary sources. BPASS WDBs in the LISA frequency contribute to almost all the signal strength that are seen by LISA, shown in Fig. D8, this could be because the ‘non-LISA’ BPASS WDBs are less massive and/or really far from us.
<figure><img src="image_25.png" /><figcaption>**Figure D8.** Same as in Fig. 15 but now including the total signal from BPASS WDBs with frequency ranging $10^{-10}-10^1$ Hz that are above the two models.</figcaption></figure>

This paper has been typeset from a `TEX/LATEX` file prepared by the author.