Neutron Star - White Dwarf Binaries: Probing Formation Pathways and Natal Kicks with LISA

Neutron star-white dwarf (NS+WD) binaries offer a unique opportunity for studying NS-specific phenomena with gravitational waves. In this paper, we employ the binary population synthesis technique to study the Galactic population of NS+WDs with the future Laser Interferometer Space Antenna (LISA). We anticipate approximately $\mathcal{O}(10^2)$ detectable NS+WDs by LISA, encompassing both circular and eccentric binaries formed via different pathways. Despite the challenge of distinguishing NS+WDs from more prevalent double white dwarfs in the LISA data (especially at frequencies below 2 mHz), we show that their eccentricity and chirp mass distributions may provide avenues to explore the NS natal kicks and common envelope evolution. Additionally, we investigate the spatial distribution of detectable NS+WDs relative to the Galactic plane and discuss prospects for identifying electromagnetic counterparts at radio wavelengths. Our results emphasise LISA's capability to detect and characterise NS+WDs and to offer insights into the properties of the underlying population. Our conclusions carry significant implications for shaping LISA data analysis strategies and future data interpretation.


INTRODUCTION
Neutron star -white dwarf (NS+WD) binaries hold significant potential in advancing our understanding of binary neutron star (NS) formation channels, natal kicks, and other NS-specific physics.In recent years, the field has experienced considerable growth particularly in the area of hydrodynamical and nuclear-hydrodynamical simulations, including the first 1-dimensional simulations of NS+WD mergers Metzger (2012), followed by significant advancements in 2dimensional simulations by Fernández & Metzger (2013) and Zenati et al. (2019), and the introduction of 3-dimensional smooth particle hydrodynamic simulations in Bobrick et al. (2022) and magnetohydrodynamic simulations in Morán-Fraile et al. (2024).These simulations have highlighted the importance of NS+WD mergers within the transient astronomy field, contributing to a deeper understanding of various transient phenomena (see also Margalit & Metzger 2016;Kaltenborn et al. 2023;Kang et al. 2023).Moreover, studies have also underscored the role of NS+WD binaries in the formation of ultra compact X-ray binaries (Tauris 2018), Thorne-Żytkow objectslike objects (Paschalidis et al. 2009), rocky planets around pulsars (Margalit & Metzger 2017), and planetary nebulae (Ablimit & Soker ★ E-mail: korol@mpa-garching.mpg.de(VK) 2024).The forthcoming Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017;Colpi et al. 2024) -as well as similar planned space-based gravitational wave (GW) observatories Tian-Qin (Luo et al. 2016;Huang et al. 2020) and Taĳi (Ruan et al. 2018) -offers a unique opportunity to discover NS+WD binaries and shed light on the underlying physics that govern their formation and evolution.In addition, in a more distant future, proposed GW missions like the Lunar Gravitational Wave Antenna (Harms et al. 2021;Branchesi et al. 2023) and the DECi-hertz Interferometer Gravitational wave Observatory (Seto et al. 2001;Arca Sedda et al. 2020) will capture GW signals from NS+WD mergers at deci-Hertz frequencies (Morán-Fraile et al. 2024).
Although neutron star -neutron star (NS+NS) systems have received more attention in the context of the LISA mission (Yu & Jeffery 2015;Kyutoku et al. 2019;Andrews et al. 2020;Lau et al. 2020;Korol & Safarzadeh 2021;Wagg et al. 2022;Seto 2022;Storck & Church 2023), it is becoming increasingly clear that NS+WD binaries may provide cleaner probes for studying NS natal kicks and other binary evolution associated phenomena (e.g.Tauris 2018;Ruiter et al. 2019;He et al. 2024).Firstly, NS+WD systems typically have a less complex evolutionary history, involving only one supernova during their formation.This simpler history makes it easier to disentangle the effects of natal kicks from other processes, allowing for a more direct assessment of the kicks' impact on the binary system.Furthermore, the relative rarity of NS+NS systems compared to NS+WD binaries means that there is a larger population of the latter available for study (e.g.Nelemans et al. 2001b;Breivik et al. 2020).By analysing a greater number of NS+WD systems, in addition to known populations of binary radio pulsars with WD companions observed and ultra-compact X-ray binaries (e.g.Rivera-Sandoval et al. 2015;Kruckow et al. 2021;Armas Padilla et al. 2023), we can obtain more robust statistics on their properties.Examining the NS+WD population's orbital characteristics, such as periods and eccentricities, can provide insights into their formation channels and the role of natal kicks (Amaro-Seoane et al. 2023).
In this work, we employ a suite of NS+WD binary population models from Toonen et al. (2018).We focus on the NS+WD parameters measurable with LISA, including the chirp mass, eccentricity, and the 3-dimensional (3D) position in the Milky Way, seeking signatures that provide insights into NS natal kicks and NS+WD formation pathways.The chirp mass -a combination of the component masses that determines the GW signal's evolution -is a crucial parameter for distinguishing NS+WD from other types of Galactic binaries accessible to LISA (mainly WD+WD and NS+NS, e.g., see Amaro-Seoane et al. 2023).In addition, we anticipate that the chirp mass can be used to distinguish between models involving different mass-transfer prescriptions.Eccentricity can also provide valuable information about the dynamical interactions and mass-transfer processes that shaped the system's orbital evolution.For example, Lau et al. (2020) have emphasised the importance of eccentricity in distinguishing various NS+NS formation channels, highlighting the crucial role of the last mass-transfer phase prior to the NS formation in determining system's characteristics.Finally, we also explore whether the observed Galactic distribution of NS+WD binaries can reveal insights into NS formation processes.Drawing from the methodology of Repetto et al. (2012Repetto et al. ( , 2017)), we also investigate if an offset of a NS+WD system from the Galactic plane can serve as a signature of peculiar velocity with respect to circular Galactic motion, expected for systems not receiving kicks.

ASSEMBLING A SYNTHETIC NS+WD POPULATION
We construct a representative present-day Galactic population of NS+WD binaries through a three-step process.Firstly, we utilise a suite of binary population synthesis models from Toonen et al. (2018) to account for the evolution from the initial main sequence (MS) stage to NS+WD formation (see Section 2.1).Each simulated catalogue represents a stellar population of approximately ∼ 10 8 M ⊙ , accounting for the fact that stars may be single or in binary systems, and may or may not have reached the main sequence turnoff within the Hubble time.These models provide binary properties at the time of NS+WD formation.Secondly, we seed these NS+WD binaries within a Milky Way potential at a rate determined by an adopted star formation history (see Section 2.2).From the NS+WD formation to the present day, we evolve the binaries' orbital parameters in accordance with GW radiation reaction (Section 2.2.1).Lastly, we assign to binaries' 3D positions in the Milky Way gravitational potential and integrate their orbits from the time of NS formation (when receiving the natal kick) until the present age of the Galaxy (Section 2.2.2).We provide a detailed description of these steps below.

From MS+MS to NS+WD
We utilise a suite of nine NS+WD population models compiled in Toonen et al. (2018) using the SeBa stellar and binary evolution module (Portegies Zwart & Verbunt 1996;Nelemans et al. 2001a;Toonen et al. 2012).Below we summarise the set of the assumptions used to generate primordial binaries and refer the reader to Toonen et al. (2018, section 3.2) for further details and discussion.Firstly, masses of primaries ( 1 ) are drawn from the initial mass function of Kroupa et al. (1993, see also Kroupa 2008) focusing on a mass range of 4 − 25 M ⊙ .We note that we consider a broader range from 0.1 to 100 M ⊙ to calculate the normalisation (i.e. the corresponding simulated stellar mass) of this population.Masses of secondaries ( 2 ) are then determined based on a flat mass ratio distribution where 0 <  2 / 1 < 1 (Raghavan et al. 2010;Duchêne & Kraus 2013).The orbital separations are drawn from a flat distribution in log(), following Abt (1983), while initial eccentricities are derived from a thermal distribution, as per Heggie (1975).Finally, we assume a constant binary fraction of 75 per cent, which is supported by observations in the considered mass range (Raghavan et al. 2010;Duchêne & Kraus 2013;Sana et al. 2014).In each model, primordial binaries and most input physics remain consistent, with variations only in the NS natal kick and common envelope prescriptions (see Sections 2.1.1 and 2.1.2).The primary goal of this study is to understand differences in the LISA-detectable NS+WD population when varying the natal kick prescription to assess whether NS natal kicks can be constrained based on the LISA data.Variations in the common envelope prescription aim to assess the uncertainty range in the LISA-detectable population, as common envelope prescriptions have been found to produce the largest differences in the population synthesis of compact binaries (e.g.Toonen et al. 2012Toonen et al. , 2014;;Korol et al. 2017;Storck & Church 2023).

NS natal kicks
It is essential to emphasise that NS formation in a binary system involves two types of kicks: (1) those due to mass loss during the supernova explosion, and (2) an additional NS-specific natal kick.The first type, referred to as 'Blaauw' kicks (Blaauw 1961), is particularly relevant in close binary systems, where interactions between the stars can strip the donor star's envelope before the supernova event.The 'Blaauw' kick magnitude can vary, depending on factors such as mass ratio, initial orbital parameters, and mass loss during the explosion.Generally, its effect is expected to be limited compared to the NS natal kick (Huang 1963;Tutukov & Yungelson 1973;Leonard et al. 1994).The formation mechanism for this additional natal kick remains unsolved (e.g.Kusenko & Segrè 1996;Scheck et al. 2006;Wongwathanarat et al. 2013;Holland-Ashford et al. 2017;Katsuda et al. 2018), but likely involves anisotropies in neutrino losses and/or mass loss in the supernova ejecta (for a review, see Janka 2012).Significant progress has been made recently in core-collapse supernova numerical simulations and resulting NS natal kicks (e.g.Coleman & Burrows 2022;Burrows et al. 2023;Janka & Kresse 2024).
NS natal kicks were discovered observationally (Lyne & Lorimer 1994) and are observed as large peculiar velocities of isolated radio pulsars (100 -1000 km s −1 ), typically at least an order of magnitude higher than the peculiar velocities of NS progenitors.These natal kicks have been investigated using proper motion measurements combined with dispersion measure and electron density models (Arzoumanian et al. 2002;Hobbs et al. 2005).Due to advancements in instrumentation -specifically, the usage of the Very Long Baseline Array with broadband phase modelling (e.g.Brisken et al. 2002) -it is now possible to measure parallax and proper motions for a large number of isolated radio pulsars (Deller et al. 2019).An assumption that peculiar velocities of young, isolated radio pulsars represent the natal kicks of NSs has several limitations as discussed by Igoshev et al. (2021) and Mandel & Igoshev (2023).In addition, there is substantial evidence that NSs formed in binaries also receive natal kicks.For instance, the number of X-ray binaries in the Small Magellanic Cloud is much smaller than expected if all NSs received only 'Blaauw' kicks (e.g., see Igoshev et al. 2021).The strength of natal kicks received by NSs formed in binaries is not well-constrained and is a subject of active research (Igoshev et al. 2021;Willcox et al. 2021;O'Doherty et al. 2023).Pfahl et al. (2002) argued that NSs in high-mass X-ray binaries with specific eccentricities and orbital periods received natal kicks below 50 km s −1 .Observational studies of orbital and spin periods for NSs in Be X-ray binaries suggest two separate populations of NSs with different natal kicks (Knigge et al. 2011).
In binary population synthesis, natal kicks are commonly modelled using the prescription developed by Fryer et al. (2012).This prescription relates the remnant mass to the initial stellar mass, and the absolute value of natal kick for core-collapse supernovae is drawn from a single Maxwellian velocity distribution with the distribution parameter  = 265 km s −1 , as estimated by Hobbs et al. (2005).The natal kicks are drawn uniformly on a sphere.The resulting orbital changes are typically calculated using equations from Brandt & Podsiadlowski (1995).
Both natal kicks and mass loss kicks contribute to a binary system's systemic velocity following a supernova explosion.We emphasise that the binary's systemic velocity represents the combined motion of both stars within the Galactic reference frame.This velocity accounts for not only the individual motion of each star but also any additional factors such as natal kicks.In the context of NS+WD binary systems, the systemic velocity is influenced by the natal kicks experienced by the NS during its formation and the system's orbital momentum.While some exceptional cases may exhibit systemic velocities of up to 600 km s −1 , typical values are expected to be around 150 km s −1 (Toonen et al. 2018), also see discussion in Section 5.2.Our binary population synthesis includes modelling of the systemic velocity.
In this study we examine four alternative supernova kick prescriptions for the NS+WD population models: • 'Blaauw' prescription: the kick is calculated as the unbalanced orbital momentum resulting from sudden mass loss (Blaauw 1961).
• 'Hobbs' prescription: the kick is drawn from a Maxwellian distribution (Hobbs et al. 2005): with the distribution parameter  = 265 km s −1 .
• 'Verbunt' prescription: Verbunt et al. (2017) determined empirically that the natal kick distribution consists of two Maxwellians with  = 0.42,  1 = 75 km s −1 , and  2 = 316 km s −1 .While this distribution represents the peculiar velocities of young isolated NS well, the physical origin of the low-and high-velocity components is still unclear.This distribution is our default choice.
Finally, kick directions are assigned isotropically.

Common envelope evolution
Observations of double compact objects in compact binaries and their mergers have prompted the concept of a phase that results in substantial orbital shrinkage (e.g.Ostriker & Davidson 1973;Paczynski 1976;van den Heuvel 1976;Meyer & Meyer-Hofmeister 1979).This phase may occur when one of the stars within a binary system evolves into a giant, and mass-transfer between the stars may happen if the binary's orbit is comparable in size to the giant star.If the material lost by the donor star (the giant) exceeds the amount the companion star (the accretor) can receive, a runaway process can develop, leading to the companion becoming fully engulfed by the giant's envelope.This unstable, runaway mass-transfer phase is referred to as the common envelope (CE) phase (for a review, see Ivanova et al. 2013;Röpke & De Marco 2023).Conceptually, the orbit of the companion star within the common envelope generates drag -either gravitational or hydrodynamical.This drag facilitates the transfer of orbital energy and angular momentum to the envelope material, thereby unbinding it from the system.As a result, the companion star spirals inward towards the core of the giant, reducing binary's orbital separation.
In binary population synthesis, the CE phase is often modelled based on energy conservation (Paczynski 1976;Webbink 1984;Livio & Soker 1988;de Kool et al. 1987;de Kool 1990).The binding energy of the envelope,  bind , is related to the loss of orbital energy Δ orb through the equation: where  d is the donor star's mass,  c is the mass of its core,  is its radius, and  is the structure parameter of its envelope that depends on the structure of the donor.The efficiency of converting orbital energy to unbind the envelope is represented by , which is often chosen based on observational calibrations (e.g.Nelemans et al. 2000;Zorotovic et al. 2010;Scherbak & Fuller 2023).As a result, the final orbital separations (and, as a consequence, the delay times between the end of the common envelope phase and the end of the gravitational wave inspiral) of the binary populations can vary depending on this choice.A higher  value corresponds to a more efficient envelope ejection.The alternative model for CE evolution, known as the -CE, focuses on the balance of angular momentum rather than that of energy (Nelemans et al. 2000).It is represented as follows: where Δ and Δ are respectively the angular momentum and mass loss.The -prescription was introduced to account for the first phase of mass-transfer during the formation of WD+WD binaries with  = 1.75 (Nelemans et al. 2000;Nelemans & Tout 2005;van der Sluys et al. 2006).
Based on these two common envelope prescriptions, we construct two families of NS+WD population models:  and .In model , the -formalism is used to determine the outcome of all CE phases.In  models, the -prescription is applied unless the binary contains a compact object or if the CE is triggered by a tidal instability (rather than dynamically unstable Roche lobe overflow, see Toonen et al. 2012).Thus, in the  models the first CE phase is typically modelled with the -CE prescription; the second CE (with a giant donor and white dwarf or neutron star companion) is typically described by the -formalism.Our default model assumes  = 2 (referred to as ) and also considers a variation with  = 0.25 (referred to as 2).Our  model, calibrated on observed WD+WDs (specifically the second mass-transfer phase, see Nelemans et al. 2000Nelemans et al. , 2001a)), contrasts with our 2 model, which is calibrated on the formation of compact WDs in binaries with M type main-sequence stars (Zorotovic et al. 2010;Toonen & Nelemans 2013;Camacho et al. 2014;Zorotovic et al. 2014).We recognise that both CE prescriptions have their own uncertainties and limitations, and they are often chosen based on the specific context of the binary system being modelled.Here we employ both to better understand their implications on the NS+WD population detectable with LISA.

From NS+WD formation to present time
The following stage in our modelling involves adjusting the properties of the simulated NS+WD binaries (such as their spatial, orbital separation, and eccentricity distributions) to resemble those of the present-day Milky Way's stellar population.This process necessitates two distinct types of integration, each carried out independently as detailed below.The first integration evolves binary orbital parameters (semi-major axis and eccentricity), which gradually decay due to GW radiation, thereby circularising the binary.The second integration involves the NS+WD orbit within the Galactic potential.The latter depends on the choice of the Milky Way's gravitational potential, while the former relies on the selection of its star formation history.

NS+WD present-day orbital parameters
In this study, we adopt a star formation history from the chemospectrophotometric model of the Milky Way by Boissier & Prantzos (1999, see also Boissier & Prantzos 2000).This model employs an 'inside-out' formation scheme for the disk and incorporates empirically and/or theoretically justified prescriptions for the star formation rate, including metallicity-dependent stellar properties.We note, however, that our binary population models assume solar metallicity, which simplifies the complex chemical enrichment history of the Milky Way.This choice is somewhat justified by the observation that, when modelling LISA-detectable Galactic WD+WD population, variations due to different CE assumptions are significantly larger than those resulting from variations in other factors, such as the initial mass function, metallicity, and binary fraction (Korol et al. 2020).We distribute NS+WD binaries in time (and space, see Section 2.2.2) according to their (spatially resolved) star formation grid, which determines the binary's age.As in Boissier & Prantzos (1999), we assume the Milky Way's age to be 13.5 Gyr.
Subsequently, we determine the present-day orbital parameters of NS+WDs-semi-major axis  and eccentricity -by accounting for GW radiation reaction from NS+WD formation until 13.5 Gyr.To do so, we numerically solve equations originally derived in Peters (1964): and where  0 is determined by the initial condition ( 0 ) =  0 , with  0 and  0 representing the binary's semi-major axis and eccentricity at NS+WD formation, and with  being the speed of light.Lastly, we exclude binaries if their formation times exceed 13.5 Gyr, if they have initiated mass-transfer (i.e., when the WD fills its Roche lobe), or if they have already merged within this time frame.We deffer the modelling of interacting NS+WD binaries to a future study.

NS+WD present-day spatial distribution
Our goal is to investigate the size and distribution of properties of NS+WDs detectable with LISA, so we can disregard the detailed structure of our Galaxy, such as the shape of its central bulge/bar region, the structure of the stellar disk including spiral arms, and the stellar halo component.We assume an analytical expression for the gravitational potential consisting of three main components: bulge, stellar disk, and dark matter halo.
To model the star formation history of the disk, as in Nelemans et al. (2004), we linearly interpolate the plane-projected star formation rate of Boissier & Prantzos (1999, SFR BP (, )), which extend up to 13.5 Gyr in time and up to 19 kpc in radius.We assume that the probability of a binary being born at a radius smaller than  follows the integrated SFR, defined as: where 0 ≤  ≤ 19 kpc is the cylindrical radius measured from the Galactic Centre.Additionally, we assume the vertical distribution of stars that follows a sech 2 function.As a result of the assumptions above, our disk stellar number density profile can be analytically described as where  d = 2.5 kpc and  d = 300 pc are respectively the characteristic scale radius and scale height of the disk (Boissier & Prantzos 1999).By integrating SFR BP (, ) up to 13.5 Gyr, we obtain a total mass of 5×10 10 M ⊙ , which is consistent with recent studies (e.g.Licquia & Newman 2015).We note that the disk scale parameters (total stellar mass, scale radius and scale heights) in our model fall within the range of values derived from observations (see Bland-Hawthorn & Gerhard 2016, for a review).We assume the distance of the Sun from the Galactic Centre to be  ⊙ = 8.5 kpc and  ⊙ = 240 km s −1 compatible with Feast & Whitelock (1997) and Reid et al. (2014).The results by GRAVITY Collaboration et al. ( 2019) are slightly different  ⊙ = 8.178±0.013stat ±0.022 syst kpc.However, we anticipate that this difference does not affect our results.In our experience this difference could result in ≈ 10 km s −1 difference which is very small in comparison to typical systematic velocities of binaries.
We model the bulge component by doubling the star formation rate in the inner 3 kpc of the Galaxy and distributing binaries as where  is the spherical distance from the Galactic Centre and  b = 0.5 kpc is the characteristic radius.Although our modeling approach is somewhat simplistic, it is based on the inside-out star formation process described by Boissier & Prantzos (1999).Consequently, in our representation of the Milky Way, we find that the median age of binaries is approximately 10 Gyr in the bulge, decreasing to around 3 Gyr at the Solar radius, and decreasing further down to about 2 Gyr at the disk's outskirts, as expected from observations (e.g.Haywood et al. 2016;Mackereth et al. 2017;Grady et al. 2020).We find the 0.0 The total rotation curve of the Milky Way is represented by the solid black line, with the contributions from individual components depicted as follows: disk (magenta dotted line), bulge (gold dash-dotted line), and dark matter halo (blue dashed line).In our simulations, the Sun's distance from the Galactic Centre,  ⊙ , and its circular velocity,  ⊙ , are set at 8.5 kpc and 240 km s −1 , respectively.integrated total mass of the bulge ( < 3 kpc) at the present time is ∼ 2.6 × 10 10 M ⊙ .
To model the density distribution of the dark matter halo we use the Navarro-Frenk-White profile (Navarro et al. 1996): where  s = 20 kpc is the scale length of the halo and  h = 0.5 × 10 7 M ⊙ kpc −3 is the halo scale density (e.g.Nesti & Salucci 2013).The total mass of the halo can be obtained by integrating Eq. ( 11) from the centre to the maximum Galactocentric radius of 100 kpc, which for our fiducial parameters yields 4.8 × 10 11 M ⊙ .
We construct a static gravitational potential based on the presentday Galactic properties following the same model as described in Eqs.(9-11) with the help of galpy package (Bovy 2015).We compute the amplitudes of individual components of the galpy potential as the following: These amplitudes have different dimensionality depending on the structure of the respective component and how exactly the aforementioned mathematical expressions are implemented in galpy 1 .These equations are summarised here to make our work reproducible.We show the resulting rotational curve and contribution of individual Galactic components to the total gravitational potential in Figure 1.We prepare the initial conditions for the Galactic orbits integration, from receiving the NS kick until the present time, as follows.We 1 See the documentation for more details https://docs.galpy.org/remind the reader that binary's birth time is determined by Eq. ( 8), while the time necessary until NS formation is provided as part of our binary population synthesis modelling.Binary's coordinates are sampled randomly following the Milky Way density profile (cf.Eqs.9-11).The initial velocities are the sum of three components: (1) the speed of the local standard of rest (LSR),  LSR , (2) a small peculiar velocity,  p , of the progenitor binary, and (3) a randomly oriented kick,  k , received by the binary due to the NS natal kick and mass loss.
The components of the binary peculiar velocity are drawn from a normal distribution,  p ∼ N (0, ), with  = 10 km/s in each direction, see e.g.Ramírez-Tannus et al. (2021) for typical velocity dispersion of massive stars.To model the components of the binary kick, we draw two angles,  ∈ [0, 2] and  ∈ [0, ], describing the orientation of the kick.These angles are drawn uniformly on a sphere i.e. we ensure the isotropic velocity distribution.The individual kick components are then computed as: Overall, the components of the total velocity (  is radial component along the axis connecting Galactic Centre and star's position and pointing away from the Galactic Centre,   is component along the Galactic rotation and   is a component along the direction perpendicular to the disk and aligned with the Galactic spin axis) can be described as the following: We integrate the orbit of a binary in the Galactic gravitational potential from NS formation until present using the Dormand-Prince integrator (Dormand & Prince 1980) as implemented in galpy using the components of Galactic potential described above.Here we assume that the gravitational potential does not evolve with time.The lengths of the integration interval is computed as 13.5 Gyr minus the time at which NS formed that we determine as the difference between the time at which the binary has been seeded in the Milky Way (cf.Section 2.2.1) and the time required to form the NS binary component (provided by SeBa, cf.Section 2.1).At the end of the integration we check if the total energy is conserved at a level better than 10 −3 of the initial energy.We do not expect that the time-dependent components of the Milky Way such as the bar and the spiral arms can significantly affect the integrated orbits because most of our binaries receive a significant systemic velocity at the moment of supernova explosion and thus they are quite separated in the phase-space from other Galactic components and thus should interact less with them.Nevertheless, these effects should be studied in more details in future publications, especially for bulge population.We also ignore the influence of the Large Magellanic Cloud -the most massive of the Milky Way's satellites -on the Galactic gravitational potential (for example, see Avner & King 1967;Conroy et al. 2021), which could possibly affect binaries' orbits in the Galactic halo.In rare cases where energy is not conserved at the required level, we disregard the orbit.These cases typically correspond to numerically challenging orbits, which may, for example, cross the Galactic Centre.The total number of disregarded orbits is usually below approximately 1 per cent of the total.We store these present-day positions of binaries in both Cartesian and Ecliptic coordinate systems and assess their detectability with LISA as detailed below.

DETECTABILITY WITH LISA
The majority of double compact object binaries within the LISA band are expected to have circular (binary) orbits.This is because of two factors.Firstly, to reach short-period orbits binaries undergo from a few to several mass-transfer phases, at least one of these phases is the CE.It is generally assumed that the CE evolution is very efficient in circularising binary's orbit.Secondly, although at a slower rate, GWs radiation also contributes to the circularisation process.However, in cases where the binary system forms in close proximity to, or even directly within, the LISA band, the binary can retain some eccentricity if the WD component forms before the NS (cf.Section 4.1.1).This residual eccentricity can be attributed to the natal kick the NS receives at the moment of the supernova explosion (cf.Section 2.1.1).
In the context of the LISA mission, to describe the GW radiation emitted by a typical quasi-monochromatic circular binary eight parameters are needed.These typically chosen to be: where A 0 is the GW amplitude, (, ) are the ecliptic longitude and latitude, respectively,  is the inclination angle,  is the polarisation angle, and  0 is an initial phase.In the circular case, the GW frequency is twice binary orbital frequency  orb while the GW amplitude is given by This is set by the source's distance  and chirp mass for component masses  1 and  2 .The chirp mass also sets the rate at which the frequency changes due to the gravitational radiation reaction: In the eccentric case, binaries emit GWs at multiple frequency harmonics; for a Keplerian orbit these are given by (with small corrections if the orbit is precessing).Each harmonic has an amplitude where the function (, ) is defined in Peters (1964).Importantly, it should be noted that in the case of a circular binary ( = 0), the function (2, 0) = 1, as expected.The rate of frequency change becomes where  () is the enhancement factor (Peters 1964) . (24)

Extracting detectable binaries from the Galactic foreground
To get an estimate of the individually detectable NS+WD binaries from the mock catalogue, we use the analysis pipeline presented in Karnesis et al. (2021).It is based on an signal-to-noise (SNR, ) evaluation using an iterative scheme for the estimate of the confusion foreground generated by Milky Way's WD+WD and NS+WD populations (see also Timpano et al. 2006;Crowder & Cornish 2007;Nissanke et al. 2012).This scheme begins with the generation of the signal measured by LISA, by computing the waveforms for each of the simulated catalogue entries, and by projecting it on the LISA arms for a given duration of the mission.After the data generation part, we begin with an iterative source subtraction process, each iteration initialised by computing an estimate of the Power Spectral Density (PSD)  n, k , for the total noise of the instrument.The index  represents the algorithm iteration number.The  n, k represents the overall noise, which includes the combined effect of overlapping unresolved GW sources and the instrumental noise, and is given by computing the running median of the data PSD.Then we process by computing the SNR   of each source  using the smoothed  n, k .The SNR  is computed as where ℎ is the true waveform template of a circular quasimonochromatic source, and the (•|•) denotes the noise weighted inner product in frequency domain, which for two real time series  and  is written as: The ' ˜' here represents the Fourier transformation, and the ' * ' the complex conjugation.The   is the one-sided noise PSD.
If   >  0 , where  0 a given chosen SNR threshold, the source is classified as resolvable, and is thus subtracted from the data.The smoothed PSD of the residual  n, k+1 is re-evaluated after iterating over the catalogue of sources, and the procedure is repeated until the algorithm converges.Convergence is reached when all sources are subtracted given the  0 threshold, or if  n, k+1 and  n, k are practically identical at all frequencies considered.At the end of this procedure, we compute the final SNR with respect to the final estimate of  n, k final , for the sources recovered.
For every binary listed in our data sets, we conduct an analysis using the Fisher information matrix (FIM) to gauge the precision of parameter extraction.It is important to note that the error assessments from the FIM analysis hold true predominantly for large SNR values ( ≳ 20, for instance Cutler 1998).Therefore, in some cases, the derived uncertainties may be underestimated.A comprehensive Bayesian parameter evaluation is essential to ascertain more realistic uncertainties for the binary parameters (e.g.Katz et al. 2022).
Finally, we should mention again that LISA is going to measure thousands of signals originating from different types of GW sources (Galactic, extra-galactic and cosmological), overlapping in time and in frequency.This demands for taking into account the correlation between the different waveforms and their corresponding parameters.Thus, a Global Fit pipeline, based on costly computational algorithms (e.g.Littenberg & Cornish 2023), is necessary.The procedure described in this section can be considered as a simulation of a computationally costly Global Fit analysis of the LISA data, and is necessary in order to keep the computational cost in acceptable levels at this exploratory stage.

Eccentricity measurement
Depending on the binary parameters, not all harmonics may reach the detectability SNR threshold, which is described in the following section.If no harmonics reach this threshold then the source is not detectable.If just one harmonic reaches this threshold then the source is detectable, but might be mistaken for a source on a circular orbit (e.g. a WD+WD binary).If two (or more) harmonics reach the threshold then they can both (all) be detected individually and, when analysed in combination, used to measure the eccentricity.Therefore, requiring that a source has at least two harmonics above the SNR threshold gives a simple, conservative criteria for classifying a binary as having a detectable (i.e.nonzero) eccentricity.
The above criterion for identifying an eccentric source is conservative in the sense that it might still be possible to measure a nonzero eccentricity even when the second (loudest) harmonic is slightly below the threshold SNR.The possibility was investigated in detail by Moore et al. (2023) who were able to classify sources with detectabile eccentricities based on the frequency of the dominant ( = 2) harmonic,  GW , and the total SNR, summed across all harmonics, .Considering NS+WD binaries, the minimum detectable eccentricity was found to be (Moore et al. 2023) This fit has been tested for 0.5 ≤  GW /mHz ≤ 10 and  ≥ 8.These results was derived using the Balrog software package that are being developed for LISA data analysis and parameter estimation for all source types: supermassive binary black hole mergers Pratten et al. (2023), WD+WD binaries (Buscicchio et al. 2019;Roebber et al. 2020;Finch et al. 2023), and stellar-mass binary black holes (Buscicchio et al. 2021;Klein et al. 2022;Bandopadhyay & Moore 2023).

Overall NS+WD population in the mHz band
Based on our modelling results (cf.Section 2) we estimate that there are between O (10 4 ) and O (10 5 ) NS+WD binaries currently emitting GWs in the LISA frequency band of (10 −4 − 10 −1 ) Hz within the Milky Way.The exact number is dependent on the assumptions made for the NS natal kick and CE prescriptions (see Table 1), as discussed in the following.We anticipate that these differences becomes more pronounced after we apply LISA's selection effects (cf.Section 4.2).Our fiducial simulation using the -CE prescription with  = 2 and the 'Verbunt' NS natal kick prescription generates 1.75 × 10 5 mHz NS+WD binaries.We highlight that this is two orders of magnitudes lower compared to the total number of WD+WD binaries in the LISA band, which is also reported in Table 1 for comparison (see also Amaro-Seoane et al. 2023).We observe a significant decrease in the number of mHz NS+WD binaries for 2 and  CE models (see also Church et al. 2006;Toonen et al. 2018).The number of NS+WDs decreases by almost 90 per cent for the 2-CE model, and by 30 per cent for the -CE model.Table 1 also indicates that the assumption regarding the NS natal kick significantly influences our results.Specifically, when comparing models using the same CE prescription and efficiency, we find that variations with the 'Blaauw' kick prescription result in the highest number of NS+WDs within the LISA band.This outcome is expected, as in this case, the kick velocities are the lowest compared to other prescriptions.Even modest natal kicks can provide enough energy for the binary to move to large distances beyond the Galactic disk or, in extreme scenarios, to completely unbind the binary from the Galaxy.Figure 2 showcases the impact of various natal kick prescriptions on the spatial distribution of NS+WD in the Milky Way.We emphasise the positions of the binaries at formation in black, and their present-day positions in colour.It is clear that even with the application of the 'Blaauw' prescription (top left), NS+WD binaries are significantly dispersed at larger distances compared to their initial spatial distribution.This effect is most pronounced for the 'Hobbs' kick prescription.We examine this in more detail in Section 4.2.3.

Main formation pathways
Figure 3 displays the present-day Galactic NS+WD population in the frequency-chirp mass parameter space.We can visually identify two distinct groups in the figure: circular binaries are shown in pale orange colour, while eccentric binaries are represented in purple.The first and largest group, which makes up 56 per cent of the population, consists of circular NS+WD binaries (pale orange) spanning a chirp mass range between 0.35 M ⊙ and 1.2 M ⊙ .The formation of circular NS+WD binaries is possible when the NS forms before the WD.In which case, mass-transfer phases following NS formation are likely to spin up ('recycle') the NS, suppress its magnetic field, and circularise the orbit.
We note that within this group there is a small number of binaries (∼0.4 per cent) with M < 0.45 M ⊙ .Such low chirp mass values indicates that the secondary is a Helium-core WD ( WD < 0.5 M ⊙ ).These binaries form via a different evolutionary channel compared to those discussed above.The initial binary typically has a large mass ratio and a wide orbit, leading the primary star to fill its Roche lobe as it evolves into a supergiant, resulting in an unstable first masstransfer phase.Ultimately, the primary star undergoes core collapse, transforming into a NS.Subsequently, when the secondary star fills its Roche lobe, a final stage of mass-transfer is initiated.This phase is typically unstable; however, in instances where it is stable, a low-mass X-ray binary is formed.After the secondary star has lost its hydrogen envelope, it becomes a Helium-core WD.We note that these systems can make up the observed population of millisecond pulsars with low-mass companions.These systems are rare in our models because binaries with less massive companions have less orbital energy and, therefore, are more likely to merge during the CE phase.A recent study by Chen et al. ( 2021) reports that when using a strong magnetic breaking (e.g.Van et al. 2019) the parameter space for formation of NSs with Helium-core WD companions becomes larger, which could potentially increase the size of this sub-population.
The second group in Fig. 3 comprises binaries with  ≠ 0. These binaries are characterised by chirp masses M > 0.8 M ⊙ .The group represents 44 per cent of the total NS+WD binary population.The formation channel involves the WD forming before the NS.When this happens, the primary star transfers mass to the secondary as it evolves off the main sequence, significantly increasing the secondary's mass.The secondary then becomes more massive than the original primary and ends its life as a NS after the primary has already evolved into a WD.The resulting binary is typically eccentric due to the natal kick received at NS formation (cf.Section 2.1.1).
We note that qualitatively the same picture holds also for 2 and  model variations.However, in 2 model variations the part of the parameter space with M < 0.6 M ⊙ becomes depleted as these binaries merge before the present time.

LISA detectable population
The identification of a Galactic GW sources based on the SNR threshold and integration time is not solely determined by the source's in-Table 1. Number of WD+WD and NS+WD binaries in the Milky Way based on various model assumptions.In particular, we consider results presented in Column A displays the total number of NS+WD in the LISA band (10 −4 − 10 −1 Hz); Column B indicates count of detectable eccentric ( ≠ 0) NS+WDs; columns C and D specify the number of NS+WDs with measurable eccentricity using two different methods (cf.Section 3.2).Specifically, column C presents results from detecting at least two harmonics of the binary with  4yr > 7, whereas column D uses the result from applying Eq. ( 27).

CE CE efficiency NS kick
In trinsic signal amplitude.Instead, it is also influenced by the specific characteristics of the entire Galactic population, which produces the unresolved Galactic confusion foreground that is added to the LISA instrumental noise, and against which the source is compared.Therefore, we combine each of our NS+WD population models with respective populations of WD+WD binaries from Korol et al. (2020), which have been generated by adopting the same CE models (i.e. and ) and the same Milky Way model.As more sources become detectable over time in our data analysis pipeline (cf.Section 3.1), they are 'removed' from this foreground, thereby reducing its amplitude (and overall noise).The methodology for assessing the detectability with LISA, as described in Section 3, enables us to account for these effects.In our analysis we set the nominal mission duration to 4 yr and the SNR threshold for detectability to  4yr > 7.
Our simulations reveal subtle yet notable variations in the number of resolved sources, depending on whether the input catalogue includes only WD+WD binaries or both WD+WD and NS+WD binaries.These variations are within a sub-percentage level.When considering WD+WD binaries only, for the  model, the number of detectable WD+WD binaries increases to 23 × 10 3 (from 22.9 reported in Table 1), representing a variation of approximately 0.52 per cent; for the  model, the count rises to 24.5×10 3 (from 24.4×10 3 ), a slight change of about 0.22 per cent.Even variations in the NS+WD population due to different assumptions on natal kicks impact the number of detectable WD+WD binaries, as highlighted in Table 1.However, we find that the foreground signal remains relatively unchanged when comparing results for WD+WD population and the population that includes both WD+WD and NS+WD binaries.Consequently, we conclude that Galactic WD+WD population constitutes the major contributor to the unresolved Galactic foreground signal, as previously assumed in the literature.
In Fig. 4, we display the characteristic strain (i.e.ℎ  = A  √︁  obs   with  obs = 4 yr) versus GW frequency for the detectable NS+WD binaries (in colour) in our fiducial model (-CE + 'Verbunt' kick).This population is compared to that of WD+WD binaries (in grey).The black solid line shows the LISA's noise level (ESA, LISA Science Requirements Document 2018), while the grey solid line represents the unresolved Galactic foreground for an assumed 4 yr observation time, which we compute self-consistently accounting for both NS+WD and WD+WD populations as described in Section 3.1.
For each NS+WD in Fig. 4, we plot the highest SNR harmonic for binaries with non-zero (true) eccentricity, while for circular binaries we plot the  = 2 harmonic (the only available one, cf.Section 3).Our fiducial model's NS+WD detections consist of 62 binaries with  = 0 and 43 with  ≠ 0; the latter ones are highlighted with a black edge contour.We also find that detectable NS+WD binaries are formedi.e., become NS+WD after the second compact object formation -at frequencies between 10 −6 Hz and 10 −4 Hz.On average, these binaries are 3.5 Gyr old: circular ones average 1.5 Gyr, while eccentric binaries tend to be around 4 Gyr old.The figure reveals that these two NS+WD sub-populations overlap in the frequency-characteristic strain parameter space, underscoring the importance of eccentricity as a distinguishing measure for inferring their NS+WD formation pathways (see Andrews et al. 2020;Lau et al. 2020;Storck & Church 2023, for similar conclusions regarding NS+NS binaries).
Across the range of analysed models, we find that in total between ∼ 80 to ∼ 350 binaries reach  4yr > 7 after 4 years of the mission's lifetime (cf.Table 1).For comparison, in Fig. 4 we also show an alternative NS+WD population model (2-CE + 'Verbunt' kick).As will be discussed subsequently, this model variation might not accurately represent the Galactic WD+WD population (cf.Section 5).
As detailed in Section 3, there are two ways we consider by which we can label a source as having a measurable eccentricity.The first method involves the source having (at least) two harmonics with SNRs above the detection threshold.The number of sources with a measurable eccentricity by this criterion are given in column C of Table 1.The second method involves the minimum detectable eccentricity criterion in Eq. ( 27) that was derived in Moore et al. (2023) based on full Bayesian parameter estimation.The number of sources by this criterion are given in column D of Table 1.To visualise the former result, we colour-code each detectable NS+WD binaries in Fig. 4 by the number of harmonics with  4yr > 7.
In our fiducial model, we identify only 5 binaries with multiple detectable harmonics.Hence, relying on the first method (column C) to categorise binaries as eccentric could result in the misclassification of the majority of genuine eccentric sources, increasing the likelihood of confusing these binaries with circular WD+WD and NS+WD binaries.Nevertheless, when applying Eq. ( 27), the majority of the eccentric NS+WD population can be recovered.This observation suggests that while analysing the Galactic population using solely circular waveforms might serve as a good initial data analysis strategy, distinguishing between WD+WD and NS+WD binaries will necessitate eccentric waveforms.
To distinguish between circular NS+WD and WD+WD binaries, it is crucial to measure the binary's chirp mass (via the measurement of  ).For the considered range, the measurement is typically possible at  ≳ 2 mHz (e.g.Korol & Safarzadeh 2021).Figure 3 illustrates that NS+WD chirp masses range between approximately ∼0.35 M ⊙ and ∼1.2 M ⊙ , while the WD+WD chirp mass distribution is expected to peak at around ∼ 0.25 − 0.4M ⊙ (depending on the model, e.g.see Korol et al. 2022), although with a long tail extending up to 1 M ⊙ .For frequencies below 2 mHz, where only GW frequency and amplitude are measured, more massive NS+WD binaries are likely to be indistinguishable from nearby, lower-mass WD+WD binaries.Yet, if the binary's eccentricity can be determined (even in the absence of a chirp mass measurement), any detected nonzero eccentricity would hint at a NS+WD binary.This is because WD+WD binaries in the LISA frequency range are expected to be circular (e.g.Nelemans et al. 2001a;Marsh et al. 2004), while NS+NS binaries are expected to be an order of magnitude less numerous (Amaro-Seoane et al. 2023).

Comparing eccentricity distributions
In this section, we examine the NS+WD eccentricity distributions across simulated models.Figure 5 displays cumulative distributions of eccentricities for: the underlying population in the LISA band (panel A); sub-populations of detectable NS+WD binaries with  4yr > 7 (panel B); and the same sub-population, but accounting for miss-classification of eccentric binaries as circular, by (re-)assigning  = 0 to binaries with only one detectable harmonic (panel C); similarly we show sub-populations with measurable eccentricities as confirmed via Eq.( 27) (panel D).
Focusing on the panel A, it is evident that both the CE (highlighted by different line styles) and the natal kick (represented by different line colours) prescriptions play a significant role in shaping the eccentricity distribution.When comparing distributions generated with the same CE efficiency (same line styles), we observe that the 'Hobbs' natal kick prescription produces more binaries with  > 0.1 than the other models.The 'Arzoumanian' and 'Verbunt' prescriptions yield similar eccentricity distributions, while for the 'Blaauw' prescription, only a small fraction of the population (< 90 per cent) has  > 0.1.We also observe that in all 2 model variations at least half of the population has  > 0.1.Panel B, which displays the same population but with selection effects applied, we observe some changes in the shape of all lines, although they seem to maintain the relative position with respect to each other (with only one exception of 2 CE prescription + 'Blaauw' kick).The median eccentricity shifts to lower values for all models.
Panel C reveals that when analysing eccentric NS+WD binaries as a collection of circular sources and then counting how many individual binary harmonics reach the detectability threshold (assuming they can be correctly associated with the same source), it becomes evident that many genuinely eccentric binaries can be mislabelled as circular.Furthermore, this criterion appears effective for binaries with an eccentricity larger than 0.1.However, distinguishing between extreme model variations should still be feasible, such as differentiating the -CE+'Blaauw' kick prescription from the 2-CE + 'Arzoumanian'/'Hobbs' natal kick prescriptions.
Lastly, when concentrating on Panel D, it becomes evident that it offers a closer representation of the sub-sample where selection effects are applied (Panel B).However, the median value appears to be biased towards marginally higher values because the eccentricity remains non-measurable for a notable fraction of the population (cf.columns B and D in Table 1).

Comparing chirp mass distributions
As previously discussed, chirp mass measurement is essential not only for characterising NS+WD binaries but also for distinguishing them from other types of Galactic double compact objects, particularly when NS+WD chirp mass distributions may overlap with those of WD+WD and/or NS+NS binaries.In Fig. 6, we compare chirp  Firstly, we observe that in all models, NS+WD chirp masses range between 0.4 M ⊙ and 1.2 M ⊙ , peaking either at ∼ 0.9 M ⊙ (for  and  model variations) or at ∼ 1.1 M ⊙ (for 2 model variations).Thus, if the chirp mass is constrained to better than 30 per cent, NS+WDs can be separated from WD+WDs, which typically have lower chirp masses.The reason for larger NS+WD chirp masses is twofold: the presence of the more massive NS component (1.2 − 1.35 M ⊙ ) and the typically more massive WD component (0.95 M ⊙ on average).
Distinguishing NS+WD from NS+NS binaries based on the chirp mass will be more challenging, especially given that the number of NS+NS binaries is expected to be in the tens, while NS+WD binaries are anticipated to be in the hundreds (Amaro-Seoane et al. 2023).For example, models from Vigna-Gómez et al. (2020) and Fig. 3 of Korol & Safarzadeh (2021) indicate that Galactic NS+NS chirp masses may peak at ∼ 1.1 M ⊙ .This implies that a chirp mass constraint of better than 10 per cent is necessary to distinguish NS+NSs from NS+WDs in our  and  models, while the majority of NS+WDs in our 2 models may be indistinguishable from NS+NSs (based on the chirp mass only).
Focusing on the left panel of Fig. 6, we observe a significant difference between chirp mass distributions in 2 models compared to / models.It is important to recall that the difference between these models stems from the varying assumptions for the CE efficiency (cf.Section 2.1.2):our 2 models are set to be less efficient in expelling the envelope, resulting in a longer CE phase, a more significant shrinkage of the binary orbit, and a higher overall number of binaries merging before the present time.Thus, if the NS+WD population can be accurately isolated from that of WD+WD, the chirp mass distribution can provide a link for to the CE efficiency parameter.
By comparing the two panels of Fig. 6, we find that the main feature of the chirp mass distributions -i.e. the location of the peakare preserved across all model variations; although sub-structures are lost due to the reduced number of binaries.Therefore, we suggest that the position of the peak of the NS+WD chirp mass distribution can be used to distinguish between model variations assuming different CE prescriptions.

Comparing 𝑍 distributions
Lastly, we investigate how the positions of NS+WD binaries in our mock Milky Way are influenced by the NS natal kicks (cf.Fig. 2).For example, Repetto et al. (2017) demonstrated that the root mean square of the height above the Galactic plane for black holes and NSs in Xray binaries can be used as a proxy to discriminate among different natal kick prescriptions.Following the same idea, we consider the positions of NS+WD binaries as measured by LISA.In Fig. 7, we compare the  distributions across the considered NS+WD models: on the left, no selection effects are applied, while on the right we show binaries with  4yr > 7.
We observe that in the top panel, both the NS natal kick prescription and the CE efficiency effects are noticeable.The CE efficiency determines the binary orbital separation at NS+WD formation and, consequently, the binary orbital energy.Depending on the magnitude of the natal kick, the binary may or may not be disrupted.Thus, both the CE and the natal kick magnitude play a role in determining the position of a binary in the Galaxy at the present day.For a fixed CE model, the 'Arzoumanian' model produces the broadest -distribution, followed by the 'Hobbs' and 'Verbunt' models; the 'Blaauw' prescription results in a significantly narrower distribution.When LISA selection effects are applied, these differences become less pronounced because only tightest, highest frequency binaries can be detected at  larger than several kpc.Still, models using the 'Blaauw' prescription remain visually narrower compared to the rest (see also Fig. 2).

DISCUSSION
In this study we assembled a suite of NS+WD populations models using the binary population synthesis technique.We assessed the size and the properties of the LISA detectable population and investigated the ways NS+WDs can be differentiated from WD+WDs in the data based on the astrophysical properties of the population: chirp mass and eccentricity distributions as well as binaries' 3D positions in the Galaxy.In particular, distinguishing between the binary types is not trivial as both WD+WD and NS+WD binaries will appear in the LISA data as quasi-monochromatic sources and will occupy a similar parameters space (cf.Fig. 4).Therefore, it will require a combination of GW measurement (such as chirp mass, eccentricity and Galactic position) measurement and electromagnetic follow up (e.g.looking for pulsar signatures).
As the primary objective of this study is to provide guidelines for the analysis and interpretation of future LISA data, rather than fitting models to the observed sample, some of the proposed model variations may not be suitable for describing both NS+WD and WD+WD populations.For instance, the assumption regarding the efficiency parameter  in CE interactions (incorporated into our binary population synthesis modelling as the product ) plays a pivotal role in determining the number of detectable binaries (cf.Table 1 and Fig. 4).While this assumption appears to align well with studies focused on post-CE binaries (Zorotovic et al. 2010(Zorotovic et al. , 2014;;Toonen & Nelemans 2013), it does not accurately replicate the observed WD+WD population (Toonen et al. 2012(Toonen et al. , 2017)).However, a recent study by Scherbak & Fuller (2023), which performs reverse modelling of WD+WD binaries' formation histories, also derived estimates of lower CE efficiency than those assumed in our fiducial model.Furthermore, it is interesting to observe that 2 set of model variations does not introduce any confusion foreground (cf.Fig. 4).An observationallydriven WD+WD population model that incorporates distributions based on observed spectroscopic samples of WD+WD candidates produces the confusion foreground of a comparable amplitude to our fiducial model (see Korol et al. 2022).
We acknowledge that our Milky Way model may be somewhat simplistic.For example, we represent the stellar Galactic disk as a single component, rather than differentiating between the thin and thick disks.Nevertheless, we have incorporated a total stellar mass and a star formation history that encompass both the thick and thin disk stellar populations.While we have not conducted model variations to specifically explore the impact of our Galactic model on the LISA detectable NS+WD population, relevant studies for NS+NS binaries by Storck & Church (2023) and for WD+WD binaries by Georgousi et al. (2023) have indicated that changes in birth positions-such as those resulting from implementing two-disk components with distinct scale parameters-do not significantly affect the resolvable frequency or eccentricity of these binaries.Furthermore, Storck & Church (2023) have found that detailed integration of Galactic orbits for NS+NS binaries is not crucial for accurately predicting the LISA-resolved NS+NS population.However, they caution that the choice of star formation history is critical when predicting LISA-visible NS+NS populations.This consideration may also be pertinent to NS+WD binaries, especially given their potentially shorter formation timescale compared to WD+WD binaries.
We also recognise that our models only captures NS+WD binaries forming in the field and do not include those binaries that may form in globular clusters.Given that globular clusters are known to host a substantial fraction of the ultra-compact X-ray binary population, it is reasonable to expect a considerable number of NS+WD binaries in these clusters.Notably, the formation rate of low-mass X-ray binaries (LMXBs), which could potentially include NS+WD binaries, is estimated to be considerably higher per unit mass in globular clusters than in the Galactic field (e.g., Heinke et al. 2003;Armas Padilla et al. 2023).Kremer et al. (2018), for example, predicts the existence of approximately 10 3 mHz NS+WD binaries in the Milky Way's globular clusters, with a few to several of these potentially detectable by LISA.

Prospects for measuring NS+WD total mass via periastron advance
Seto (2001) pointed out that for eccentric binaries, a measurement of the periastron advance -analogous to the perihelion advance of Mercury -may be detectable in binaries in the LISA data.In the absence of periastron precession, the GW radiation is emitted at exactly multiples of the orbital frequency (Peters 1964).If the binary is precessing, the spacing of the individual harmonics will be shifted by integer multiples of the periastron precession frequency From the above equation it is clear  PP depends on the total mass  of the binary, while   depends on the chirp mass (cf.Eq. 23); Moore et al. (2023) illustrated how these two measurements can be combined to derive the individual component masses.It seems likely that for the majority of binaries where the eccentricity can be measured, it will be possible to obtain the measurement of the component masses by this technique.While this will depend on the binary's SNR and frequency, we can speculate that results reported in column D of Table 1 are representative of the number of individual component mass measurements.In the example investigated by Moore et al.
(2023, see their fig.4), the fractional uncertainty on the primary mass (90 per cent uncertainty) is ∼33 per cent.
For completeness, we highlight that Tauris (2018) discussed a method for measuring NS masses in NS+WD binaries leveraging a well-established correlation between the orbital period and the mass of a Helium-core WD in LMXB systems.Only those post-LMXB NS+WD binaries with less than about 9-hour orbital periods, which can coalesce within a Hubble time to become visible LISA sources, are observed to have WD masses within a remarkably narrow range ( WD = 0.162 ± 0.005 M ⊙ ).This fact in combination with LISA's chirp mass measurement facilitates an accurate determination of NS masses with a precision within estimated precision of ∼ 4 per cent (see their figure 4).

Systemic velocities of millisecond radio pulsars
Recently O'Doherty et al. ( 2023) performed analysis of systemic velocities for binaries with NS including millisecond radio pulsars.They integrated orbits of systems back in time and recorded their velocities at the expected moment of compact object formation.These velocities are typically smaller than the NS natal kick and ranges from tens to ≈ 250 km s −1 .Here we perform simple comparison between systemic velocities of NS+WD binaries with small eccentricities ( < 0.1) and their results.For this comparison, we select all binaries which satisfy only single aforementioned criteria i.e. systems plotted are neither selected on the base of the LISA detectability nor their radio detectability as radio pulsars.We show the result of our analysis in Figure 8.The curve obtained by O'Doherty et al. (2023) is between  and 2 models with 'Verbunt' natal kicks.The  model with Hobbs natal kick provides similar results as 2 'Verbunt' natal kicks.The small systemic velocities ( syst < 60 km s −1 ) follows 2+'Verbunt' (similarly  + 'Hobbs') cumulative distribution, but later on the curve suggested by O' Doherty et al. (2023) diverges from these curves strongly and it follows +'Verbunt' model quite closely for  syst > 120 km s −1 .
It is interesting to note that estimates of O' Doherty et al. (2023) fall between two curves produced with the same natal kick prescription but with different assumptions about the CE evolution.This suggests that values of  for millisecond pulsar formation should be between 0.25 and 2. We can also speculate that the expected LISA population would be more similar to our simulations with the 'Verbunt' natal kick prescription.

Follow up radio observations
Some of NS+WD binaries discovered by LISA could host millisecond or normal radio pulsars.It is much more probable that follow-up  radio observations discover a new millisecond radio pulsar than a normal radio pulsar because radio emission from normal pulsar is beamed much stronger than in the case of millisecond radio pulsar (e.g.Kramer et al. 1998).Having discovered NS+WD LISA source, it is possible to design an optimal strategy for radio follow-up.LISA's sky localisation for Galactic stellar-remnant binaries largely depends on the binary's SNR and frequency, and is expect to be in the range of several to few deg 2 at  GW > 2−3 mHz (e.g.Finch et al. 2023;Moore et al. 2023;Colpi et al. 2024).A typical pulsar survey covers 1.5 • in a single pointing like the Parkes multi-beam pulsar survey (Manchester et al. 2001).Thus, if NS+WD system is subsequently discovered in a dedicated pulsar search, radio observations will enormously improve the sky localisation because pulsar coordinates are typically known with accuracy exceeding 1 arcsec.Pulsar discovery will also help to constrain the orbital parameters such as the mass function and eccentricity via the pulsar timing.
The normal radio pulsar could be formed if the NS is formed after WD formation via the reversed channel.Additionally, the system should not be too old because normal pulsars stop operating typically after 0.1-1 Gyrs.Thus, highly eccentric systems with NSs can be the primary targets for dedicated young pulsar searches.The expected pulsar orbital periods will be ranging from a few tens of milliseconds to a couple of seconds, see Igoshev et al. (2022) for recent estimates on initial NS spin periods.Only a fraction of about 10 per cent of all young pulsars in NS+WD can be detected (Tauris & Manchester 1998) because in multiple cases the beam will miss the Earth due to the geometric orientation of the NS.Thus, given the number counts in Table 1 and age distribution we could expect to discover ≈ 2 − 6 young pulsars in a dedicated radio search if an efficient technique to deal with orbital accelerations is implemented.The lower estimate corresponds to the pulsar lifetime of 0.1 Gyr while upper one corresponds to lifetime of 1 Gyr.Additional constrains might arise from the limited sky coverage of radio surveys.Thus, the Parkes Multibeam Pulsar Survey (Manchester et al. 2001) covered only the Galactic longitude || < 5 • .So, these radio pulsars could be missed in radio surveys if the relevant part of the sky is not covered.
The millisecond radio pulsar is typically formed as a result of stable mass-transfer from secondary to NS.Thus, we expect that millisecond pulsar candidates have nearly circular orbit with eccentricity well below 0.1.The mass is transferred together with angular momentum which spins up the NS and it typically has a period shorter than ≈30 milliseconds.The millisecond radio pulsars could operate for many Gyr and their average radio profiles are much wider than that is for normal radio pulsars (Kramer et al. 1998).The profile width reaches 100 • .This opening angle could correspond to up to 50 per cent possible detection.Thus, many more (up to ≈ 50) are expected to be discovered in radio follow-ups among NS+WD GW emitters with small eccentricity.The caveat is that radio signal from millisecond pulsars located in regions with high free electron density (e.g. the Galactic Centre) could be smeared and leave these pulsars undetectable.

SUMMARY AND CONCLUSIONS
In this study, we examine the detectability of Galactic NS+WD binaries with LISA, estimating a population size of O (102 ) detectable systems.According to our fiducial binary population synthesis model about a half of the NS+WD population in the LISA band are circular.These are formed through the reversed channel (i.e.WD forms first, followed by the NS; see also He et al. 2024 on how this result may change when changing assumptions in the underlying binary evolution model).The remaining population consists of eccentric NS+WD binaries that originate from the direct channel, where the NS forms first, followed by the WD.In the current study, we have not undertaken a full reconstruction of the population from the simulated LISA data.Instead, our investigation has concentrated on exploring how NS+WD binaries can be differentiated from other types of Galactic binaries within the LISA dataset, as well as on identifying key observable features that are essential for aiding data inference.Given the high completeness of LISA's sample up to frequencies of 2-3 mHz (Amaro-Seoane et al. 2023), we anticipate that it will be possible to recover the underlying true formation rates and population properties from LISA data, contingent upon the identification of NS-WD binaries within the dataset.
The identification of NS+WD binaries within the LISA data may be challenging, especially at frequencies below 2 mHz, where they might be misclassified as WD+WD binaries.In particular, this misclassification risk is heightened for circular NS+WD binaries, which can be easily confused with WD+WDs when the chirp mass measurement is absent (typically at frequencies lower than a few mHz).However, eccentric NS+WD binaries can potentially be pinpointed through the eccentricity measurement.We explored two methods to identify sources with measurable eccentricity in our catalogues.The first requires a source to present at least two harmonics above the detection threshold, mimicking the analysis of the Galactic population using only circular waveforms.In contrast, the second leverages the minimum detectable eccentricity criterion, as derived in Moore et al. (2023) using eccentric waveforms for source parameters recovery.Relying exclusively on the first method carries a risk of misidentifying the majority of genuine eccentric sources.However, the second method, which employs the minimum detectable eccentricity criterion, proves more adept at distinguishing eccentric sources.While initial analyses -particularly within the LISA Global fit framework (e.g.Littenberg & Cornish 2023) -might rely on circular waveforms, a more clear differentiation between WD+WD and NS+WD binaries calls for the inclusion of eccentric waveforms (for a detailed analysis, refer to Moore et al. 2023).We aim to further investigate source confusion through a more realistic Global fit type of analysis in a future study.
We have also demonstrated that the NS+WD chirp mass distribution is sensitive to the CE efficiency parameter by comparing  and 2 sub-sets of models.In particular, we propose the idea that the position of the peak of the chirp mass distribution can be used to dis-tinguish between these model variations.Additionally, our analysis underscores the sensitivity of the eccentricity distribution to the NS natal kick prescription, emphasising the potential to distinguish (at least) between extreme kick prescriptions based on the subset overall population accessible with LISA.However, learning about the NS natal kicks based on NS+WD positions above the Galactic plane -as it has been done for high-mass X-ray binaries in Repetto et al. (2017) -might be challenging.This is due to the constraints on position measurements with LISA, a marked difference when compared to high-mass X-ray binaries, where measurements are typically more precise (deg 2 versus arcsec 2 ).
Lastly, we discussed the prospects for identifying electromagnetic counterparts of NS+WD binaries in the radio regime, which could offer additional insight into the formation channels and characteristics of these systems.Our preliminary calculations suggest the possibility of detecting a few young pulsars with WD companions and several tens of millisecond pulsars through targeted radio follow-up.
In conclusion, our exploration of the detectability of Galactic NS+WD binaries with the upcoming LISA mission has highlighted the potential of unlocking of NS-specific science with the LISA data.We have pinpointed challenges in classifying NS+WD binaries accurately, owing to the abundance of WD+WD binaries and the potential overlap with NS+NS binaries.Our conclusions, together with those in the companion paper by Moore et al. (2023), carry significant implications for defining LISA data analysis strategies and data interpretation.

Figure 2 .
Figure 2. Comparison of present-day NS+WD positions in  −  Galactocentric coordinates under varying NS natal kick prescriptions, depicted in different colours.The initial (MS+MS) positions are indicated in black for reference.

Figure 3 .
Figure 3. Displayed is our fiducial model of the present-day Galactic NS+WD population ( -CE + 'Verbunt' kick), represented in the frequency -chirp mass parameter space (bottom left), with projected frequency and chirp mass distributions per bin presented in the top left and bottom right panels respectively.Binaries are categorised by colour: circular ( = 0, pale orange) and eccentric ( ≠ 0, purple).The top right panel showcases the distribution of eccentricities.The dashed grey line in the top left panel shows the expected ∝  −11/3 distribution in frequency, with fewer sources at higher frequencies because of the accelerating, or chirping, inspiral described in Eq. (23).

Figure 4 .
Figure 4. LISA-detectable Galactic NS+WD (in colour) and WD+WD (in grey) populations in the GW frequency-characteristic strain parameter space.The top panel presents our fiducial model ( -CE with  = 2 and the 'Verbunt' natal kick), while the bottom panel displays an alternative model ( 2-CE with  = 0.25 and the 'Verbunt' natal kick).We only plot binaries that have SNRs of  4yr > 7. NS+WD binaries with  ≠ 0 are highlighted with a black contour.For these eccentric binaries, individual harmonics must meet the SNR criterion.For clarity, only the highest SNR harmonic for each binary is shown, with the colour indicating the number of detectable harmonics specific to that binary.The black solid line showcases LISA's sensitivity curve (ESA, LISA Science Requirements Document 2018).The grey solid line illustrates the unresolved Galactic foreground, calculated self-consistently from the two populations (NS+WD and WD+WD binaries) for a 4-year mission duration.

Figure 5 .
Figure 5. Cumulative eccentricity distributions for NS+WD populations, referenced to Table 1 columns A, B, C and D. Panel A (top left): underlying Galactic NS+WD population within the LISA band.Panel B (top right): sub-populations with  4yr > 7. Panel C (bottom left): sub-populations with more than one detectable harmonic.Panel D (bottom right): sub-populations with measurable eccentricities, as assessed using Eq.(27).Different line styles (solid and dashed) represent models assuming different CE efficiency parameter, while colour variations indicate different natal kick prescriptions, consistent with Fig. 2. Distributions in panels B, C, and D are normalised to the total number of detectable binaries in their respective models.

Figure 6 .
Figure 6.Comparison of the chirp mass density distributions across the suite of considered models.Left panel: no selection effects; right panel: including selection effects.

Figure 7 .
Figure 7.Comparison of NS+WD distribution in  coordinate across the suite of considered models.Left panel: no selection effects; right panel: including selection effects.

Figure 8 .
Figure 8. Distribution of NS+WD systemic velocities for binaries with small eccentricities  < 0.1.