Discovering neutron stars with LISA via measurements of orbital eccentricity in Galactic binaries

LISA will detect $\sim \! 10^4$ Galactic binaries, the majority being double white dwarfs. However, approximately $\sim \! 1 \textrm{--} 5 \%$ of these systems will contain neutron stars which, if they can be correctly identified, will provide new opportunities for studying binary evolution pathways involving mass reversal and supernovae as well as being promising targets for multi-messenger observations. Eccentricity, expected from neutron star natal kicks, will be a key identifying signature for binaries containing a neutron star. Eccentric binaries radiate at widely-spaced frequency harmonics that must first be identified as originating from a single source and then analysed coherently. A multi-harmonic heterodyning approach for this type of data analysis is used to perform Bayesian parameter estimation on a range of simulated eccentric LISA signals. This is used to: (i) investigate LISA's ability to measure orbital eccentricity and to quantify the minimum detectable eccentricity; (ii) demonstrate how eccentricity and periastron precession help to break the mass degeneracy allowing the individual component masses to be inferred, potentially confirming the presence of a neutron star; (iii) investigate the possibility of source misidentification when the individual harmonics of an eccentric binary masquerade as separate circular binaries; and (iv) investigate the possibility of source reclassification, where parameter estimation results of multiple circular analyses are combined in postprocessing to quickly infer the parameters of an eccentric source. The broader implications of this for the ongoing design of the LISA global fit are also discussed.


INTRODUCTION
The Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) is a planned space-based gravitational-wave (GW) observatory being developed for launch in the mid 2030s.LISA will observe the source-rich mHz GW spectrum, which is populated by ∼ 10 8 shortperiod (∼ 1 hour) stellar-mass compact-object binaries in our galaxy (Lipunov et al. 1987;Hils et al. 1990;Nelemans et al. 2001b) and its satellites (Roebber et al. 2020;Korol et al. 2021).These weakfield GW sources generate quasi-monochromatic signals and will be the most numerous LISA sources.The majority will be double white-dwarf (WD+WD) binaries; in addition to a few tens of known verification systems (see, e.g., Kupfer et al. 2023;Finch et al. 2023) LISA is expected to discover a few tens of thousands of individually detectable WD+WD binaries (Nelemans et al. 2004;Nissanke et al. 2012;Kremer et al. 2017;Breivik et al. 2020;Li et al. 2020;Korol et al. 2020Korol et al. , 2022) ) and to be able to probe the stochastic GW signal ★ E-mail: cmoore@star.sr.bham.ac.uk generated by the entire Galactic population of several million binaries (Hils et al. 1990;Bender & Hils 1997;Farmer & Phinney 2003;Ruiter et al. 2010;Georgousi et al. 2023).
LISA will also detect Galactic quasi-monochromatic binaries including a neutron star (NS) (Korol 2023) and/or black hole (Wagg et al. 2022) component, albeit in far smaller numbers.Of particular interest here are the several hundred neutron star -white dwarf (NS+WD) binaries (Breivik et al. 2020;Korol et al. 2024) that are now expected to be detectable by LISA (this is significantly more than early estimates, Nelemans et al. 2001a).Formation of these sources is expected from isolated binaries evolving in the field with a significant fraction of binaries undergoing mass reversal (Toonen et al. 2018;Tutukov & Yungelson 1993;Tauris & Sennels 2000;Davies et al. 2002;Bobrick et al. 2017;Zenati et al. 2019) during a period of stable mass transfer which results in the WD forming before the NS.Then, when the NS forms, the accompanying supernovae kick puts the binary onto an eccentric orbit.This picture is supported by observations of several binary-pulsar systems with WDs on eccen- The harmonic frequency structure for an illustrative set of six eccentric quasi-monochromatic binaries.Each binary is either low frequency (  GW = 0.5 mHz, blue shades) or high frequency (  GW = 6 mHz, green shades) and either low ( = 0.08, down-pointing triangles), medium ( = 0.2, circles), or high (0.25, up-pointing triangles) eccentricity.Harmonics arise at half-integer multiples of the frequency of the dominant  = 2 harmonic.Also shown is the LISA noise curve, with and without the expected contribution from the Galactic confusion noise for a 4 year mission.For each harmonic, the height of the marker relative to the noise curve indicates the SNR in that harmonic.The harmonics are extremely narrow band; the inset plot shows the structure of one harmonic in one LISA TDI channel.The spikes clearly visible in the inset are spaced at frequency intervals of 1 year −1 and are due to the orbital modulation of the instrument response.All sources were simulated for a 4-year LISA mission, with  GW = 5.99 × 10 −19 Hz s −1 ,  PP = 3.91 nHz,  = /4, and in the direction of the Galactic centre (ecliptic coordinates of [, ] = [4.66,−0.098] radians).The total SNR (summed across all harmonics) of each source was chosen to be  = 10.tric orbits (Lyne & McKenna 1989;van Kerkwĳk & Kulkarni 1999;Kaspi et al. 2000;Brown et al. 2001;Ng et al. 2018).
Observing NS+WD binaries with LISA will help in understanding the mass transfer and supernovae kicks in their evolution as well as providing new targets for multi-messenger astronomy (Korol 2023).An important question is if it will be possible to tell if a binary contains a NS from its GW signal alone?This is challenging because quasi-monochromatic GW signals contain very little information.For exactly monochromatic binaries, the binary masses cannot be inferred at all, and even if the binary is chirping (and assuming this is due to GW emission), only the chirp mass combination can be inferred (Korol & Safarzadeh 2021).Orbital eccentricity might be an important smoking gun; WD+WD binaries are expected to have been efficiently circularised by several episodes of mass transfer, so measuring a non-zero eccentricity strongly suggests a NS+WD source, or another rare type of Galactic binary (e.g., involving a black hole) or a WD+WD binary from an unusual evolutionary pathway.Furthermore, for binaries on eccentric orbits, the rate of advance of the periapsis in general relativity (GR) depends on the total mass of the binary; therefore, if the periastron precession (also known as apsidal precession) frequency can be measured from the GW signal then this allows the total mass of the binary to be inferred (Seto 2001, see also Eq. 12).For binaries where the chirp and periapsis precession can both be measured, it will be possible to infer the total mass and the chirp mass separately, breaking the degeneracy and allowing the individual component masses to be determined.
The effects of orbital eccentricity imprint themselves on the GW signals of quasi-monochromatic binaries in several ways (Peters & Mathews 1963;Peters 1964;Damour & Deruelle 1985;Moreno-Garrido et al. 1994, 1995;Seto 2001;Damour et al. 2001Damour et al. , 2005;;Mishra et al. 2015).The most striking feature is the presence of narrow-band harmonics at multiples of the orbital period, , with frequencies   ∼ / (see, e.g., Fig. 1).The relativistic effect of periastron precession (cf. the orbit of Mercury) causes the orientation of the orbital ellipse to advance over time and precess with a frequency  PP ≪  =2 .This periastron precession modulates the GW signal, further shifting the frequency of the harmonics by multiples of  PP (an expression for the frequency of the  th harmonic is given in Eq. 22).Finally, orbital eccentricity also enhances the overall flux of GWs, leading to a quicker inspiral and the frequency of each harmonic increasing faster with time.
Narrow-band signals can be efficiently analysed using a heterodyning approach focusing on the small number of relevant frequency bins.This paper uses a multi-harmonic heterodyning algorithm for the efficient analysis of multiple narrow bands widely separated in frequency.A new problem of source confusion also arises when analysing eccentric signals.The multiple, narrow frequency harmonics generated by an eccentric source can be mistaken for separate circular sources.It is possible to accurately match an eccentric GW source with a series of circular templates (one for each harmonic), thereby mistaking one source for several.This has implications for LISA data analysis where the numerous sources in the millihertz band generate GW signals that are overlapping in both time a frequency space and must all be analysed together in a process called a global fit.Different parts of the global fit algorithm designed to target different source types (and which may need to be run asynchronously) must be tolerant to mistaken source identifications made by other parts of the algorithm.Different parts of the algorithm will need to share information to find and correct such mistakes, e.g. by reclassifying multiple circular source candidates as a single eccentric source.
The structure of this paper is as follows: The basic physics of eccentric Galactic binaries and their GW signals is described in Sec. 2. The data analysis methods used are described in Sec. 3. The three main results of this paper are presented in Sec.4: Sec.4.1 quantifies LISA's minimum measurable eccentricity for quasi-monochromatic Galactic binaries as a function of source parameters; Sec.4.2 verifies that for certain systems, eccentricity allows the component masses to be inferred, potentially confirming the presence of a NS; Sec.4.3 considers the problem of source misidentification; and Sec.4.4 considers the problem of source reclassification.Finally, Sec. 5 concludes by discussing some of the broader implications of this work for ongoing work on designing the LISA global fit.

GW PHYSICS OF ECCENTRIC BINARIES
Quasi-monochromatic Galactic binaries are weak field sources of GWs.Therefore, the GW signals can be computed (in the transversetraceless gauge) using the quadrupole formula (Misner et al. 1973).This has been investigated for a binary on a Keplerian elliptical orbit by several authors (see, e.g., Moreno-Garrido et al. 1995;Seto 2001).We employ the following expressions for the two polarizations: (1) This signal contains discrete harmonics at frequencies   , with corrections coming from the periastron precession.Here, A  () is the amplitude of each harmonic,  is the inclination,  is the phase, Φ is the argument of periastron, and  is the eccentricity.
Typically, the number of harmonics  included should be chosen depending on the eccentricity and the frequency in order to capture a given fraction of the squared signal-to-noise ratio (SNR; see Eq. 19).However, to facilitate comparisons between sources with different eccentricities, most of the calculations in this paper were performed with a fixed value  = 12.For the small to moderate eccentricities considered in this paper, this is a conservative choice as there is negligible power at larger .The effect of doing the analysis with different numbers of harmonics is investigated in Appendix C.
Instead of the orbital period , it will be more convenient for us to work with the initial frequency of what is (usually, at least when  ≲ 0.29; see Fig. 2) the dominant harmonic,  GW =  =2 ( = 0) = 2/.(The time  = 0 is defined as the start of LISA science operations.)A binary radiating GWs does not stay on a fixed Keplerian orbit and to account for this the frequency is allowed to drift, or chirp, linearly in time according to The argument of periastron is also allowed to advance according to where  PP is the periastron precession frequency.Note that we use different conventions than those in the references, which explains the different form of our polarizations.For example, when comparing with Moreno-Garrido et al. (1995) -Garrido et al. 1995).We also choose to give the frequency of each harmonic in terms of the periastron-to-periastron frequency,  peri (such that  GW = 2  peri ), whereas the references give frequencies relative to the orbital frequency  orb (that is, the frequency of sweeping out 2 radians).These are related via  orb =  peri +  PP .
Our GW polarization amplitudes, A  , are implemented via an expansion correct up to tenth order in eccentricity.While this approach is obviously not suitable for eccentricities close to unity, it has been found to work extremely well for the small to moderate eccentricities  ≲ 0.6 expected in WD+NS binaries.To make the connection to Moreno-Garrido et al. (1995), the harmonic amplitudes are the expansions of where A is an overall amplitude common to each harmonic (see Eq. 14).For low values of eccentricity, we can also write where (, ) is a function proportional to the power in the  th harmonic (notation from Peters & Mathews 1963).When expanded about  = 0, Eqs. 4 and 5 are equal up to and including the ( + 1) th power in  (for  ≥ 2).
The LISA response to GWs must also be modelled.This depends on the sky position (ecliptic coordinates  and ) and orientation (via the polarisation angle ) of the source.The orbital motion of LISA modulates the response at a frequency of 1 year −1 and introduces Doppler broadening into each frequency mode.The LISA phase measurements are processed into three time series containing the GW information through a process known as time-delay interferometry (TDI; Tinto & Dhurandhar 2021) which are designed to suppresses the large-amplitude laser noise below the level of the other noise sources.The model for the LISA response is described in Klein (2021).Schematically, , ,  = LISA_Response(ℎ + , ℎ × ; , , ), where ,  and  are the noise-orthogonal TDI output channels.The log-likelihood is expressed in terms of these time series (see Eq. 16).
Eqs. 1 to 6 constitute our GW waveform model for eccentric quasimonochromatic binaries.There are 11 independent signal parameters: {A, ,  GW ,  GW ,  PP , , , Φ 0 , , , }.When working with the cosine of the inclination and the sine of the ecliptic latitude, inference was performed directly on these parameters with uniform priors.We chose to allow the parameter inference to sample values of the parameters that would be nonphysical if the binary was simply two point masses moving under the influence of gravity alone (for example, we are free to use priors that allow negative  GW or  PP ).This is done to allow the analysis to detect binaries that may be interacting, e.g.via mass transfer or tidal forces.The conversion to the physical binary parameters described below is performed, if required, in post processing.It should be remembered that these conversions will give biased results (or may fail entirely) if the binary is interacting.
Consider a binary with masses  1 and  2 and semi-major axis .The orbital period is related to  and the total mass  =  1 +  2 by Kepler's third law; in terms of  GW , this is As the binary inspirals the frequency increases, or chirps, with time GW 1 mHz 11/3  () × 0.18 nHz year −1 . (10) Here, an overdot denotes differentiation with respect to time, M  = ( 1  2 ) 3/5 /( 1 +  2 ) 1/5 is the binary chirp mass, and is the eccentricity enhancement (see Fig. 2; Peters & Mathews 1963).At leading post-Newtonian order, the argument of periastron advances with time.The periastron precession frequency is The frequency derivative depends on M  while the periastron advance depends on .Therefore, if both  GW and  PP can be measured, it is possible to determine  1 and  2 independently.Finally, if the source is at a distance , the GW strain amplitude is given by

LISA DATA ANALYSIS WITH MULTIPLE HARMONICS
The most important ingredient for Bayesian GW data analysis is the likelihood function.Under the common simplifying assumptions that the noise in LISA has a known power spectral density (PSD) and is additive, Gaussian and independent in each of the three noiseorthogonal TDI channels, the likelihood is given by Here,  is the data, ℎ() is the waveform model as a function of the 11 source parameters, , and angled brackets denote the signal inner product in each of the three TDI channels,  ∈ {, ,  }.
The factors not shown in the proportional relation in Eq. 16 are unimportant constants in the sense that they do not depend on .
The inner products are defined as integrals (or sums for discretely sampled data) over frequency.However, the long LISA mission and accompanying fine frequency resolution (Δ  = 1/ obs ) coupled with quasi-monochromatic nature of these sources mean that it would be extremely inefficient to sum over all frequencies.Instead, the signal inner product in Eq. 16 is split into a sum of integrals targeting narrow frequency ranges around each harmonic; where A subscript  denotes the time series in the relevant TDI channel and a star denotes complex conjugation.The integration limits (  ,low ,  ,high ) are chosen to target each frequency harmonic based on the source priors, and the bandwidths |  ,high −  ,low | are chosen to be wide enough to capture all the signal power in a given harmonic after accounting for the source chirp and Doppler broadening caused by the motion of the LISA constellation, while avoiding spectral leakage between different harmonics; this can occur when negative frequency components of one harmonic get heterodyned onto another harmonic.For example,  ,high is determined by taking the upper end of the prior ranges for  GW and  PP combined using Eq.22, plus a term depending on the upper end of the prior range for  GW calculated using Eq.23, plus a term corresponding to the maximum possible Doppler shift if the source happens to lie in the ecliptic plane.Each of the narrow-band sums in Eq. 18 can then be efficiently evaluated by using heterodyning to shift the frequency range of interest close to zero, low-pass filtering the data and down sampling to target only the small number of frequency bins in the range  ,low to  ,high (with several extra padding bins included for safety).
An important quantity in GW data analysis is the (optimum) SNR.This largely determines if a source is detectable and any parameter estimation uncertainties.The SNR is defined as The total source SNR is the quadrature sum of the SNRs in the individual harmonics; i.e.  2 =  =1  2  , where The inner products in Eq. 16 were evaluated using a simple analytic model for the noise PSD based on the latest LISA science requirements (SciRD).This model includes estimates of the contributions from both the instrumental noise (following Babak et al. 2021) and the Galactic confusion noise (using Eq. 4 of Babak et al. 2017 To a good approximation, the minimum eccentricity of a quasi-monochromatic binary that LISA can measure depends only on the frequency of the source and its SNR.These curves were constructed using a 4-year LISA mission, with source properties corresponding to a 1.4  ⊙ -1.4  ⊙ binary in the direction of the Galactic centre, at an inclination  = /4; however, as verified in Appendix A, the results are not sensitive to these choices.The inset shows a sequence of eccentricity posteriors for a 3 mHz source at an SNR of 13, with four different injected eccentricities.The significance at which each eccentricity posterior is peaked away from zero (the median divided by the standard deviation of the distribution) is indicated along the vertical lines (which are placed at the injected value of the eccentricity).This shows visually what is meant by the threshold criterion of requiring the median over standard deviation to exceed 3 before claiming a eccentricity detectable.These results are confirmed using instead a threshold on the evidence ratio (or Bayes' factor) in Appendix B.
the population of eccentric sources (not present in Babak et al.) does not contribute to the confusion foreground.The analysis was performed using the Balrog software package that is being developed for LISA data analysis and parameter estimation for all source types.This includes supermassive binary black hole mergers Pratten et al. (2023), double WDs (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).
Parameter estimation was performed using the nessai implementation (Williams 2021;Williams et al. 2021) of the nested sampling algorithm (Skilling 2004(Skilling , 2006)).Sampling was performed directly on the 11 parameters described in Sec. 2 using flat priors.

The Minimum Detectable Eccentricity
This section aims to establish the minimum eccentricity,  min , that a quasi-monochromatic binary must have in order for LISA to be able to measure it.In other words, the smallest eccentricity that allows us to confidently conclude a binary is not circular.It is expected that  min will depend strongly on  GW and the SNR, , as these directly control the strength of the various harmonics relative to the LISA noise curve.The influence of the other parameters is expected to be smaller and/or tend to average out after several orbits of the LISA constellation; this is verified in Appendix A.
The analysis framework described in Sec. 3 was used to perform a series of zero-noise injections with different values of  GW ,  (which was controlled by varying A), and .Sources were simulated at 13 values of SNR,  ∈ {7, 8, . . ., 19}, and 16 values of frequency,  GW /mHz ∈ 0.5, 0.75, . . ., 2.75, 3, 4, 5, 6, 8, 10}.For each pair of  GW and , a series of analyses were performed at different  to find  min .All other source parameters were held constant at the values described in the caption of Fig. 3.The significance of the eccentricity measurement was quantified via the ratio of the median to the standard deviation, denoted , of the 1-dimensional marginal posterior on  (see inset plot in Fig. 3).For simplicity, the eccentricity is considered detectable if this ratio exceeds a fixed threshold value of 3.These results were also checked using the Bayes' factor between an analysis using circular and eccentric waveform models with similar results (see Appendix B).
The minimum detectable eccentricity is plotted in Fig. 3.For convenience, a fitting formula for these results is also provided; This fit has been tested for 0.5 ≤  GW /mHz ≤ 10 and  ≥ 8.For large , the fit scales as  min ∝  −1 , so Eq.21 can be extrapolated to high SNRs.However, this formula should not be extrapolated to lower SNRs, where the  −1 dependence breaks down and sources can become hard to detect at all.Indeed, at  = 8, the maximum residual and error of Eq.21 is 0.015 and 14% respectively, whereas for all  ≥ 9 this improves to 0.0055 and 11% respectively.Outside of the tested frequency range the LISA sensitivity is expected to decrease quite rapidly.
The minimum detectable eccentricities in Fig. 3,  min ∼ 0.1, are fairly large in comparison to other types of sources in LISA.The difficulty in measuring eccentricity is a consequence of the small amount of information contained in the quasi-monochromatic GW signals.This contrasts strongly with LISA's ability to measure much smaller eccentricities in other, broadband signals.For example, when The top panel shows how the measurements of the chirp mass (coming from the chirp rate  GW ) and the total mass (coming from the periastron precession  PP ) combine to break the mass degeneracy and allow both  1 and  2 to be determined separately.Vertical dot-dashed lines indicate the approximate upper and lower ends of the WD and NS mass distributions in a particular population synthesis model (described in the text) and show that from the mass measurements (together with that the binary is eccentric;  = 0.1 ± 0.005) we could confidently conclude that this binary contains a NS from the LISA observation of its GW signal alone.The bottom row shows selected posteriors for similar binaries with smaller eccentricities of  = 0.05 (red) and  = 0.01 (purple), with SNRs 29.3 and 27.2 respectively.(When  = 0.01, the measured eccentricity is consistent with zero and it is no longer possible to determine the component masses.) observing supermassive binary black hole mergers, eccentricities as small as  ≲ 10 −2 are expected to be measurable (see, for example, Garg et al. 2023).Or, for chirping stellar mass black holes near the upper end of the LISA sensitive band, eccentricities as small as  ≲ 10 −3 are expected to be measurable (see, for example, Klein et al. 2022).Or, for extreme-mass-ratio inspirals, eccentricities as small as  ≲ 10 −4 are expected to be measurable (see, for example, Barack & Cutler 2004;Babak et al. 2017).
The optimum frequency where LISA's sensitivity to eccentricity in quasi-monochromatic binaries is best is around ∼ 2 mHz.This is because this places the third harmonic (which is usually the second loudest) at a frequency of ∼ 1.5 × 2 = 3 mHz near the minimum of the LISA noise bucket (see Fig. 1).Measuring the first subdominant harmonic (usually  = 3) confirms that the source is eccentric.
The actual number of NS+WD systems detectable by LISA and the fraction of those that will have a measurable eccentricity by this criterion can be assessed using binary population synthesis.This has been done by Korol et al. (2024) where it is shown that (using their fiducial model) there are 105 LISA-detectable NS+WD binaries of which 25 have a eccentricity that is measurable by LISA (as assessed using Eq.21).These numbers are rather uncertain; the range of population models considered give a range 75-357 for the total number of detectable NS+WD systems and 6-68 for those with measurable eccentricities as assessed using Eq.21 (Korol et al. 2024, Table 1).

Determination of the Component Masses
It will generally not be possible to determine the component masses of quasi-monochromatic binaries observed by LISA from their GW signals alone.Even if the binary has a detectable chirp,  GW , that we are confident is due solely to the emission of GWs, only the chirp mass combination M  can be determined (see Eq. 9).
High-eccentricity quasi-monochromatic binaries are an exception to this.This section considers the possibility of using the measured GW chirp and the periastron precession frequency to determine the individual component masses of the binary.This is possible because  GW depends on the chirp mass (see Eq. 9), whereas  PP depends on the total mass (see Eq. 12), and measuring both breaks the  1 ,  2 parameter degeneracy.Measuring the individual component masses and determining that at least one is above ≳ 1.35  ⊙ would be the most direct way of confirming the presence of a NS in a quasimonochromatic binary with LISA alone.
Here it is assumed that the binary contains two point masses interacting only gravitationally with no external or environmental effects.In this situation, the measured signal parameters {A,  GW ,  GW ,  PP } can be related to the physical binary parameters { 1 ,  2 , , } as described in Sec. 2. It should be emphasised that if these assumptions fail to hold, biased parameter estimates will be obtained.
There are several possible cases to be considered: (i) If the binary is exactly monochromatic and only  GW is measured, then only the Kepler parameter combination / 3 can be inferred (Eq.7).This reveals little about the nature of the binary.(ii) If  GW can also be measured, then the chirp mass M  can be inferred (Eq.9); see Korol & Safarzadeh (2021).However, knowing M  alone doesn't unambiguously determine if the binary contains a NS.(iii) If  and  PP (but not  GW ) can also be measured, then the total mass  can be inferred (Eq.12); see Seto (2001).The eccentric orbit suggests the binary evolved through a supernova, which is suggestive of a NS.However, knowing  alone doesn't unambiguously determine if the binary contains a NS.(iv) If  GW ,  GW ,  and  PP can all be measured, then it is possible to infer the orbital semi-major axis  and the individual component masses  1 and  2 separately.Together with the fact that the binary is eccentric, this is the strongest possible evidence from GWs alone that the binary contains a NS.
Fig. 4 shows parameter estimation results for a binary in case (iv).Because not all combinations of the signal parameters correspond to a physical, non-interacting binary; when converting to the binary parameters, additional prior cuts have to be imposed (see discussion in Sec. 2).This is achieved in practice by discarding posterior samples with nonphysical parameters combinations; for the  = 0.1, 0.05, 0.01 binaries in Fig. 4, the fractions of discarded samples were 0.28, 0.43 and 0.47 respectively.The posteriors on the physical parameters in Fig. 4 were plotted with the remaining physically valid samples.If, in a real analysis, the majority of samples are nonphysical, this may indicate that the source is not an isolated, non-interacting binary.
The top panel of Fig. 4 is designed to show how the simultaneous measurements of the chirp mass (from the measurement of ) and the total mass (from the measurement of  PP ) break the mass degeneracy.The primary mass is determined to be greater than 1.3  ⊙ with 90% confidence; strongly suggesting that this is a NS.
Fig. 4 also includes vertical lines separating regions of the primary mass range where we expect either WDs or NSs.These boundaries are determined using the fiducial NS+WD mock population in Korol et al. (2024) which was assembled using binary population synthesis.In this model the WD component mass probability density peaks ∼ 0.9  ⊙ with a tail extending up to the Chandrasekhar limit (∼ 1.4  ⊙ ), although such high masses are significantly less likely and the vast majority of primary objects are below 1  ⊙ .The minimum NS mass in NS+WD binaries is ∼ 1.22  ⊙ .Comparing these boundaries to observations of WD+WD systems poses challenges, particularly in constraining the mass of the primary star (the more massive of the pair).This task can be complicated by the secondary star, which is often brighter and, consequently, presents difficulties in characterising the primary's properties.Nevertheless, we can compare these boundaries to the mass distribution of single WDs; in spectroscopic samples, WDs with masses exceeding 1  ⊙ are highly improbable and likely to be a result of a merger (Kepler et al. 2015;Temmink et al. 2020;Kilic et al. 2023).The distribution of NS masses in binaries has a mean at ∼ 1.28  ⊙ with a dispersion of 0.24  ⊙ (Özel et al. 2012, see also Farrow et al. 2019).
For the other source parameters not shown in Fig. 4, orbital eccentricity has only a minor impact on our ability to measure them compared to the case of circular quasi-monochromatic binary.The parameters that show the largest effect are the sky position.For the binary in Fig. 4 with  = 0.1, LISA can localise the source to a 90% sky area of Ω 90 = 7.9 deg 2 whereas for the same binary with  = 0.01 this area increases to Ω 90 = 16.5 deg 2 .Much of this difference can be attributed to the increased SNR of the eccentric source (34.8 compared to 27.2), the expected scaling being Ω 90 ∝  −2 .The remainder of the difference is likely due to the presence of some signal power at higher frequencies for the eccentric binary, where LISA's angular resolution is known to be better (Cutler 1998).

Misidentification of Harmonics as Circular Sources
Each harmonic of an eccentric quasi-monochromatic binary appears in the LISA data as a narrow-band feature in the frequency spectrum.
It is possible (in fact, very easy) to mistake these features as being produced by a number of separate circular quasi-monochromatic binaries, where the number of such circular sources is simply the number of harmonics with an SNR above the detection threshold (approximately  ≥ 6).This sort of source confusion has implications for the design of the LISA global fit (see Sec. 5).
Using the waveforms described in Sec. 2, an eccentric source that generates  harmonics above the SNR threshold can be (mis)modelled perfectly by  circular sources, meaning the overlap between the circular sources and the eccentric harmonics is unity.The (fictitious) circular templates that give this perfect overlap all have the same sky position and inclination as the (true) eccentric source.The initial frequency of the circular template matching the  th harmonic of the eccentric source is given by Similarly, the chirp rate of the circular template matching the  th harmonic of the eccentric source is given by Finally, the amplitude of the circular template matching the  th harmonic of the eccentric source is given by Eqs. 4 or 5.In general, the  = 2 harmonic has the largest amplitude (at least for moderate eccentricities,  ≲ 0.29) and this is a decreasing function of eccentricity, while the amplitudes of the other modes are smaller and are increasing functions of the eccentricity.This behaviour can also be seen clearly in Fig. 1.
In order to investigate the possibility of confusing eccentric harmonics for circular sources another, zero-noise injection was performed using a loud eccentric source with  = 0.3 and  = 21.7.Parameter estimation was first performed on this signal using the (correct) eccentric waveform model.Then, parameter estimation was performed again using an (incorrect) circular waveform model targeting only the data in a narrow frequency range around each harmonic.A subset of the parameter estimation results are shown in Fig. 5.For this system, the SNRs in the first five harmonics were  =1 = 0.273,  =2 = 14.6,  =3 = 13.7, =4 = 7.42, and  =5 = 3.44.It was found that the circular analyses could only recover the  = 2, 3, and 4 harmonics, since these are the harmonics above the detectability threshold of ∼ 6.
As expected, the eccentric analyses correctly recover all the injected source parameters.The circular analyses recover parameter values that can be related back to the parameters of the eccentric injection.For example, the frequency (bottom panel) of the  th harmonic is a function of the frequency of the eccentric injection,  GW,ecc , its periastron precession frequency,  PP , and the harmonic index,  (see Eq. 22).The sky positions (top right) of all the (fictitious) circular sources are consistent with that of the eccentric injection but with larger uncertainties corresponding to the reduced SNR when only a single harmonic is being analysed.Likewise, the inclinations (top left) of the circular sources are consistent with the eccentric injection but with larger uncertainties.The amplitudes (top left) of the circular sources are lower than that of the eccentric injection, again corresponding to the reduced SNR in a single harmonic and given by Eq. 4.
In a global fit scenario, the simpler (and more common) circular sources will likely be targeted first, and so this misidentification of eccentric harmonics as separate circular sources is likely to occur (this is discussed further in Sec. 5).Inclination, ι / rad Posteriors illustrating the possibility of confusing eccentric harmonics for separate, circular sources.The eccentric source used for the injection had parameters  GW,ecc = 3.0 mHz,  PP = 310 nHz, A = 1.3 × 10 −23 ,  = /4,  = 0.3 and was in the direction of the Galactic centre.The total SNR of the injection was  = 21.7.The individual harmonics of this eccentric injection can be mistaken for separate, circular binaries.Results are shown for four independent parameter estimation analyses: recovery with an eccentric model (ecc, the truth), recovery with a circular model focusing on frequencies around ∼ 3.0 mHz (c 2 , the  = 2 harmonic of the injected source), recovery with a circular model focusing on frequencies around ∼ 4.5 mHz (c 3 ), and recovery with a circular model focusing on frequencies around ∼ 6 mHz (c 4 ).Additional circular analysis were performed targeting the  = 1 and  = 5 harmonics; however, nothing was recovered because these harmonics were too quiet.In all 2-dimensional plots, pairs of iso-probability contours enclose 50% and 90% of the posterior probability.We emphasise that the results in this figure were obtained with zero-noise injections.This was done for simplicity.If realistic noise were to be included the effect would be to shift the posteriors of the frequency for each the individual circular binaries up or down by an amount consistent with the posterior width.This will make the process of identifying these circular sources as harmonics of a single eccentric source slightly harder than it appears here; a detailed analysis of this problem involving a large number of noisy injections is beyond the scope of this paper.Eccentric analysis (ecc) Inferred from c 2 and c 3 analyses Figure 6.Posteriors on selected parameters for the same eccentric source analysed in Fig. 5.The pale dashed line shows the results from the eccentric analysis (using the same waveform model for the parameter estimation analysis as was used for the injection).The darker solid line shows the rough results that can be inferred by combining the results of two circular parameter estimation analyses targeting the  = 2 and  = 3 harmonics (the c 2 and c 3 analyses described in Fig. 5).The three parameters shown are the eccentricity , the periastron precession frequency  PP and the strain amplitude, A. The spacing of the recovered frequencies from the two circular analyses (the fact that the frequency of the c 3 source is not exactly 3/2 times the c 2 source) encodes  PP and the ratio of the recovered amplitudes from the two circular analyses encodes .In all 2-dimensional plots, pairs of iso-probability contours enclose 50% and 90% of the posterior probability.

Reclassification of Circular Sources as an Eccentric Source
As mentioned in the previous section, the first stage of a global fit is likely to involve searching for all circular sources that can be found in the data, as these are the most numerous and simplest LISA sources.In this scenario, a search for eccentric sources will only be performed at a later stage.As part of this search, it will be possible to identify pairs of circular source candidates that are likely two harmonics of a single, eccentric source (e.g. they have overlapping sky and inclination posteriors and with frequencies that are closely related by a factor of 3/2).Using the existing circular parameter estimates for this pair of sources, it is possible to get a crude estimate for the eccentric source parameters without any additional parameter estimation; the frequency spacing gives you the periastron precession frequency, and the amplitudes give you the eccentricity.This can be used to verify if two circular candidates are in fact harmonics of a single eccentric source (by seeing if the inferred  PP and  are astrophysically reasonable, for example).In a global fit the multiple circular waveforms will want to be swapped out for a single eccentric waveform, and this procedure could also be used to efficiently initialise the parameters of the eccentric waveform.This section demonstrates the idea of estimating eccentric parameters from the circular parameter estimation for the eccentric source in Fig. 5.For simplicity we use only the circular analyses of the  = 2 and  = 3 harmonics.
Let  GW,c n denote the value of the frequency of a circular template being used to analyse the  th harmonic of an eccentric source.The relationship between the circular frequencies and the properties of the eccentric source is given in Eq. 22.This was used with  = 2 and  = 3 to infer the value of  PP for the eccentric source.The strain amplitudes A c  of the circular templates can be related to the eccentricity.The expressions for the individual amplitudes of the  = 2 and  = 3 harmonics are given by Eq. 4. These expressions were inverted numerically to find A and  as functions of A c 2 and A c 3 .The resulting posteriors are shown in Fig. 6.
As can be seen from Fig. 6, the posteriors on the eccentric parameters inferred from the circular analyses are broader (i.e. less informative) than those from the eccentric analysis.This is expected because the circular analyses are not analysing all of the available signal coherently, and this results in a loss of SNR.The SNRs in the individual harmonics for this source are given in Table 1.

DISCUSSION AND IMPLICATIONS FOR THE DESIGN OF THE LISA GLOBAL FIT
A few percent of the quasi-monochromatic binaries that LISA will observe may be eccentric.A non-zero eccentricity suggests that the source may contain a NS and merit followup analyses.If both an increasing frequency and advancing periastron can be measured in a non-interacting binary then the individual component masses can be determined, helping to confirm the presence of a NS.This has been demonstrated for several simulated sources and LISA's ability to measure eccentricity has been quantified across parameter space.
Using the results for the LISA minimum measurable eccentricity in combination with binary population synthesis, Korol et al. (2024) expect between 6-68 of the 75-357 detectable WD+NS LISA sources to have a measurable eccentricity.The widely-separated frequency harmonics of eccentric quasimonochromatic binaries are unusual among GW sources, which generally consist of overlapping broadband modes.A multi-harmonic heterodyning approach that allows all these harmonics to be efficiently analysed coherently has been used as part of our Bayesian parameter inference.The benefits of analysing all harmonics together (including harmonics that contain a small fraction of the total SNR), as opposed analysing harmonics individually as if they were circular sources has been demonstrated (see Appendix C).
Ultimately, both eccentric and circular sources will have to be analysed together as part of the LISA global fit.This consists of the simultaneous detection and characterisation of the numerous GW sources, instrumental features, glitches, and noise sources in the LISA data.The details of this complex process are yet to be finalised (for an early prototype showing what this might look like, see Littenberg & Cornish 2023).However, it is clear that a hierarchical and iterative approach will be necessary, with the simplest sources being provisionally identified and analysed first before the rest are tackled.This process will iterated many times before it eventually converges.
The first stage of the global fit will almost certainly target the easily-resolvable circular (i.e.single-harmonic) Galactic binaries as these are both the simplest and most numerous LISA source type.As has been shown here, this first stage may incorrectly identify a harmonic of a loud eccentric source as a circular source.(If the eccentric source is sufficiently loud and eccentric, then several harmonics might be incorrectly identified as separate circular sources.)It will be necessary for the later stages of the global fit to correct these misidentifications, either by looking for circular source candidates in the still-evolving catalogue that are likely harmonics of a single source (e.g., with frequencies that differ by a factor of approximately /2 and that have consistent sky positions, frequency derivatives and inclination angles) and removing/replacing them, or by regularly performing eccentric versus circular model comparison on every candidate in the continually-evolving catalogue.This highlights a general feature required by any LISA global fit; while the algorithm as a whole is still running, the inputs and outputs of its different components are only provisional and must be subject to revision and modification by other parts of the algorithm.
It is also interesting to note that, what with the large number of circular sources in the LISA band, it is possible that multiple (actually circular) sources could by chance have a frequency spacing that makes them appear as a false-eccentric source.Multiple sources with overlapping sky and inclination posteriors is not surprising, since the uncertainty on these is large compared to the separation of sources on the sky.Frequency is the best measured and will be the best distinguisher.However, note that there is some freedom in the spacing introduced by the periastron precession factor, which will not be known beforehand.The combined phase-polarization measurement could also be useful, but only if more than two modes are detected (the phase-polarization combination for each harmonic will be regularly spaced, but that spacing is not known beforehand so this is not useful when only two harmonics are measured).And note that this will be a bigger problem earlier in the mission, when less SNR has accumulated and posteriors are wider.
Eccentric quasi-monochromatic binaries are a possible LISA source type that has hitherto received little attention in the literature.While, like circular Galactic binaries, these sources are relatively simple to model they nevertheless pose a new analysis challenge due to the necessity of analysing multiple, widely separated frequency harmonics and the possibility of confusion with the more numerous circular sources.This paper has shown how these problems may be tackled and how, by identifying eccentric sources, LISA may be used as a new tool for the discovery of NS+WD Galactic binaries.Finally, the dependence of the minimum detectable eccentricity on the mission duration was investigated.The fiducial source was reanalysed with simulated mission durations of 1, 2, . . ., 10 years, see Fig. A1.The SNR increases with time, leading to improved eccentricity measurements.For mission durations longer than a couple of years, the improvement in the eccentricity measurement has the expected scaling of  −1 , indicating that there is no other effect on the eccentricity measurement besides that of the increased SNR.
Taken together, these results show that LISA's ability to measure eccentricity in quasi-monochromatic binaries depends only very weakly on all the source and mission parameters except the source frequency,  GW , and the SNR, .This justifies the use of the fitting formula in Eq.21 that is a function of just these two parameters.

APPENDIX B: USING THE BAYES' FACTOR TO QUANTIFY THE MEASURABILITY OF ECCENTRICITY
The results for the minimum detectable eccentricity,  min (  GW , ), in Fig. 3 were obtained by requiring that the 1-dimensional marginalised posterior on the eccentricity parameter was peaked sufficiently away from zero.This was quantified using the ratio of the median to the standard deviation of the posterior distribution on  and this ratio was required to be greater than three in order to claim the eccentricity as measurable.
While this threshold for eccentricity detection has the advantage of being simple to implement and simple to interpret visually (see inset in Fig. 3) it is also somewhat arbitrary and difficult to relate to other statistical measures of detection confidence.Here we reproduce similar results for  min using a threshold on Bayes' factor.The Bayes' factor is defined as the ratio of the Bayesian evidence between an analysis using an eccentric waveform model and an analysis using a circular model ( = 0).An equal prior odds ratio was used between these two models.The Bayes' factor is calculated in two ways: using the ratio of the two evidences calculated via nested sampling and (exploiting the fact that the circular model is nested in the eccentric model) using the Savage-Dickey density ratio on the parameter .
Using the Bayes' factor as a statistic for detecting the eccentricity, we now claim a measurement of the eccentricity if the Bayes' factor exceeds a certain predetermined threshold, .The results are shown in Fig. B1.This analysis gives an independent check on the results in Fig. 3 and allows us to establish that our previous threshold of three on the ratio of median to the standard deviation corresponds to a threshold on the Bayes' factor of  = 2.3.

APPENDIX C: THE EFFECT OF INCLUDING SUBDOMINANT HARMONICS
Most of the parameter estimation calculations performed so far in this paper have either been circular analyses targeting a single harmonic, or eccentric analyses targeting all  = 12 harmonics.This appendix investigates eccentric analyses that analyse only a smaller number of harmonics with  in the range 2 to 12.
The benefit of analysing a small number of harmonics is increased computation efficiency.The disadvantage is the loss of SNR when only analysing a limited number of harmonics.Although this is mitigated by the fact that the SNR is concentrated in a small number of harmonics for the small to moderated eccentricities expected for Galactic WD+NS binaries (see Table 1).
The eccentric source used in Fig. 5 was reanalysed using a different number of total harmonics.The results are shown in Fig. C1.As expected, the effect of including a particular harmonic scale with the SNR in that harmonic.However, even harmonics that contain small SNR (significantly below the threshold of ∼ 6 that would be required for them to be identified as an individual source in a circular source search) still have a noticeable effect on the posterior.This paper has been typeset from a T E X/L A T E X file prepared by the author.the analytic fit in Eq.21 to the data in Fig. 3 obtained using a threshold on the median over standard deviation of the posterior on .Circular (star-shaped) markers show the results obtained using a threshold on the Bayes' factor calculated using the nested sampling (Savage-Dickey) method.These were obtained from a sequence of 5 runs with increasing eccentricity (at fixed values of all other parameters), computing the Bayes' factor for each run, linearly interpolating the Bayes' factor as a function of eccentricity, and determining the value of  min where the Bayes' factor exceeds a threshold .The top panel shows the Bayes' factors as a function of eccentricity for an example sequence of runs at  GW = 3 mHz.The threshold on the Bayes' factor of  = 2.3 was determined by minimising the root-mean-square (r.m.s.) difference between the circular points and the dashed lines, the resulting r.m.s.difference is Δ min = 0.006.If instead the r.m.s.difference between the Savage-Dickey results (star-shaped points) and the dashed lines is minimised, a threshold of  = 1.9 was obtained, with an r.m.s.difference of Δ min = 0.004.Overall, the Bayes' factor results agree well with the results in Fig. 3, justifying the use of the detection threshold used.
Figure1.The harmonic frequency structure for an illustrative set of six eccentric quasi-monochromatic binaries.Each binary is either low frequency (  GW = 0.5 mHz, blue shades) or high frequency (  GW = 6 mHz, green shades) and either low ( = 0.08, down-pointing triangles), medium ( = 0.2, circles), or high (0.25, up-pointing triangles) eccentricity.Harmonics arise at half-integer multiples of the frequency of the dominant  = 2 harmonic.Also shown is the LISA noise curve, with and without the expected contribution from the Galactic confusion noise for a 4 year mission.For each harmonic, the height of the marker relative to the noise curve indicates the SNR in that harmonic.The harmonics are extremely narrow band; the inset plot shows the structure of one harmonic in one LISA TDI channel.The spikes clearly visible in the inset are spaced at frequency intervals of 1 year −1 and are due to the orbital modulation of the instrument response.All sources were simulated for a 4-year LISA mission, with  GW = 5.99 × 10 −19 Hz s −1 ,  PP = 3.91 nHz,  = /4, and in the direction of the Galactic centre (ecliptic coordinates of [, ] = [4.66,−0.098] radians).The total SNR (summed across all harmonics) of each source was chosen to be  = 10.

Figure 4 .
Figure 4. LISA parameter estimation results for a highly eccentric WD+NS binary (green contours).The simulated binary parameters were:  1 = 1.4  ⊙ and  2 = 0.9  ⊙ ,  = 16.7 minutes (meaning  GW = 2 mHz),  = 0.1,  = /3,  = 10 kpc, and in the direction of the Galactic centre.The simulated binary is detached, meaning  GW and  PP take the GR values given by Eqs. 9 and 12 respectively.These parameter choices give a total source SNR of  = 34.8 after 4 years of observation.Joint 2-dimensional posterior distributions are shown on the parameters  1 ,  2 ,  and .(The convention  1 ≥  2 excludes the grey region of parameter space.)The top panel shows how the measurements of the chirp mass (coming from the chirp rate  GW ) and the total mass (coming from the periastron precession  PP ) combine to break the mass degeneracy and allow both  1 and  2 to be determined separately.Vertical dot-dashed lines indicate the approximate upper and lower ends of the WD and NS mass distributions in a particular population synthesis model (described in the text) and show that from the mass measurements (together with that the binary is eccentric;  = 0.1 ± 0.005) we could confidently conclude that this binary contains a NS from the LISA observation of its GW signal alone.The bottom row shows selected posteriors for similar binaries with smaller eccentricities of  = 0.05 (red) and  = 0.01 (purple), with SNRs 29.3 and 27.2 respectively.(When  = 0.01, the measured eccentricity is consistent with zero and it is no longer possible to determine the component masses.)

Figure B1 .
Figure B1.The bottom panel shows the minimum detectable eccentricity as a function of SNR, for a subset of the GW frequencies tested.Dashed lines show the analytic fit in Eq.21 to the data in Fig.3obtained using a threshold on the median over standard deviation of the posterior on .Circular (star-shaped) markers show the results obtained using a threshold on the Bayes' factor calculated using the nested sampling (Savage-Dickey) method.These were obtained from a sequence of 5 runs with increasing eccentricity (at fixed values of all other parameters), computing the Bayes' factor for each run, linearly interpolating the Bayes' factor as a function of eccentricity, and determining the value of  min where the Bayes' factor exceeds a threshold .The top panel shows the Bayes' factors as a function of eccentricity for an example sequence of runs at  GW = 3 mHz.The threshold on the Bayes' factor of  = 2.3 was determined by minimising the root-mean-square (r.m.s.) difference between the circular points and the dashed lines, the resulting r.m.s.difference is Δ min = 0.006.If instead the r.m.s.difference between the Savage-Dickey results (star-shaped points) and the dashed lines is minimised, a threshold of  = 1.9 was obtained, with an r.m.s.difference of Δ min = 0.004.Overall, the Bayes' factor results agree well with the results in Fig.3, justifying the use of the detection threshold used.
, we only keep dominant (  +   )/2 terms (where   and   are functions of  involving the Bessel functions of the first kind   ; see Eqs.A15 and A16 of Moreno the The improvement in eccentricity measurement with increasing mission duration,  obs .Top panel: posteriors on the eccentricity.Middle panel: the posteriors on the periastron precession frequency.Blue (orange) horizontal lines indicate the injected parameter values and the median and 50% credible intervals are indicated in each violin plot.The eccentricity can be confidently measured to be greater than zero for all  obs > 4 years and the posterior widths on both  and  PP scale as  −1 for longer mission durations.Bottom panel: the total signal SNR.For a 1 year mission the SNR is below 6 (gray shaded region) and the source cannot be confidently detected (this is why no posteriors on  or  PP are plotted for  obs = 1 year).