How can LISA probe a population of GW190425-like binary neutron stars in the Milky Way?

The nature of GW190425, a presumed binary neutron star (BNS) merger detected by the LIGO/Virgo Scientiﬁc Collaboration (LVC) with a total mass of 3 . 4 + 0 . 3 − 0 . 1 M (cid:2) , remains a mystery. With such a large total mass, GW190425 stands at ﬁve standard deviations away from the total mass distribution of Galactic BNSs of 2.66 ± 0.12 M (cid:2) . LVC suggested that this system could be a BNS formed from a fast-merging channel rendering its non-detection at radio wavelengths due to selection effects. BNSs with orbital periods less than a few hours – progenitors of LIGO/Virgo mergers – are prime target candidates for the future Laser Interferometer Space Antenna ( LISA ). If GW190425-like binaries exist in the Milky Way, LISA will detect them within the volume of our Galaxy and will measure their chirp masses to better than 10 per cent for those binaries with gravitational wave frequencies larger than 2 mHz. This work explores how we can probe a population of Galactic GW190425-like BNSs with LISA and investigate their origin. We assume that the Milky Way’s BNS population consists of two distinct subpopulations: a fraction w 1 that follows the observed Galactic BNS chirp mass distribution and w 2 that resembles chirp mass of GW190425. We show that LISA ’s accuracy on recovering the fraction of GW190425-like binaries depends on the BNS merger rate. For the merger rates reported in the literature, 21 − 212 Myr − 1 , the error on the recovered fractions varies between ∼ 30 and 5 per cent.


I N T RO D U C T I O N
GW190425 is a compact object merger with a total mass of 3.4 +0.3 −0.1 M that was recently detected by the LIGO/Virgo Scientific Collaboration (LVC; Abbott et al. 2020). If GW190425 is a binary neutron star (BNS) merger, its total mass is inconsistent with the observed Galactic BNS population that shows a narrow range in the total mass of ≈2.66 ± 0.12 M (Farrow, Zhu & Thrane 2019). LVC suggests that GW190425 might belong to a class of BNSs born from a fast-merging channel. In this hypothesis, massive BNSs like GW190425 merge on time-scales of 10-100 Myr due to being formed with short orbital periods through unstable case-BB mass transfer or with high eccentricities through large natal kicks (e.g. Tauris et al. 2017;Galaudage et al. 2020;Romero-Shaw et al. 2020). Short lifetimes and severe Doppler smearing, which affects shortperiod systems, make these binaries invisible to radio telescopes and thus hard to find in our Galaxy (Cameron et al. 2018;Pol et al. 2020).
However, the comparable merger rate of GW190425 with R GW190425 = 460 +1050 −390 yr −1 Gpc −3 and that of GW170817 with R GW170817 = 760 +1740 −650 yr −1 Gpc −3 derived by LVC can be challenging to account for through a fast-merging channel hypothesis for two reasons. First, massive neutron stars are expected to form from more massive progenitor stars; thus, if adopting the initial mass function of Salpeter (1955), the expected merger rate of GW190425-like systems should result in a lower merger rate for GW190425 compared E-mail: korol@star.sr.bham.ac.uk to GW170817 (for a more detailed discussion, see Safarzadeh, Ramirez-Ruiz & Berger 2020). Secondly, if fast-merging BNSs form through unstable case-BB mass transfer (e.g. Dewi & Pols 2003;Ivanova et al. 2003), they should constitute 10 per cent of the total BNS number according to a suite of simulations studied in Safarzadeh et al. (2019). In addition, a large fraction of the BNSs formed through fast-merging channels could challenge the BNS origin for the rprocess enrichment of the ultrafaint dwarf galaxies (e.g. Komiya et al. 2014;Matteucci et al. 2014;Safarzadeh & Scannapieco 2017;Safarzadeh et al. 2019).
Many binary population synthesis studies investigated the observed properties of Galactic BNSs (for an overview, see Tauris et al. 2017). Although these studies are able to reproduce the broad characteristics of the BNS population, they seem to have difficulties in matching the mass distribution of the radio population and simultaneously account for the unusually high mass of GW190425 (e.g. Tauris et al. 2017;Vigna-Gómez et al. 2018;Kruckow 2020;see also Fig. 3). Binary population synthesis models using probabilistic remnant mass and kick prescriptions form heavier BNSs; still half of their population has masses larger than those of known Galactic BNSs (Mandel et al. 2021). Thus, it is not yet clear whether more GW190425-like BNSs could be associated with the fast-merging channel or whether they require different explanations. For example, massive BNSs may not be fast merging at all; instead, the lack of radio detections could be attributed to a weak neutron star magnetic dipole moment. If pulsars are born with either a very strong or extremely weak magnetic dipole moment, they will migrate into the graveyard of pulsars and become invisible to radio telescopes ).
Regardless of the origin, if massive short-period GW190425-like binaries exist in the Milky Way, the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) will be able to detect them and -for sufficiently short orbital periods -accurately measure their chirp masses. In this work, we test whether LISA can identify massive BNSs as a distinct Galactic subpopulation. We assume that two subpopulations of BNSs reside in the Milky Way: one that follows the chirp mass distribution of known Galactic BNSs, and another that resembles GW190425. The question at hand is that given the expected merger rate of BNSs in the Milky Way, can LISA detect short-period massive BNSs and recover their true fraction?
The idea that LISA can study the existence of fast-merging channel BNS system has been introduced in other works. Recently, Andrews et al. (2020) showed that even by adopting a pessimistic Milky Way merger rate for Galactic BNS systems, LISA would detect tens of such BNSs within 4 yr of observations. Lau et al. (2020) arrived at similar conclusions suggesting that with the achievable high accuracy on the chirp mass and eccentricity, one can constrain BNS formation scenarios such as the natal kicks imparted to neutron stars at birth. In this work, we quantify how many events it would take for LISA to set bounds on the fraction of GW190425-like systems in the Milky Way. It is important to mention that in this work we do not make any assumptions regarding BNS GW frequencies. If GW190425like systems come from the fast-merging channel, it is possible that they would appear at higher frequencies in the LISA band compared to the bulk of the Milky Way population. Alternatively, they could form at similar frequencies as the currently known Galactic BNSs but with high eccentricities. Thus, we anticipate that to associate massive BNSs with the fast-merging channel, additional measurements of the frequency and eccentricity distributions would be required.
The structure of this work is as follows. In Section 2, we analyse LISA's capability to determine the fractional error on the chirp mass of binaries as a function of the gravitational wave (GW) frequency. In Section 3, we set up a multivariate Gaussian distribution for the binaries detectable by LISA, one following the chirp mass distribution of the Galactic binaries and the other resembling GW190425. We study how LISA can recover the relative fraction of the two subpopulations. In Section 4, we discuss our results and conclude.

Chirp mass measurement
To start, we have to quantify the accuracy within which LISA can measure the chirp mass. The chirp mass of a binary system is defined as where m 1 and m 2 are the primary and secondary neutron star masses, respectively. It determines how fast the binary's GW frequency f = 2/P (with P being binary orbital period) changes during the in-spiral phasė f = 96 5 π 8/3 GM c 3 where G and c are, respectively, the gravitational constant and the speed of light. Therefore, chirp mass can be derived from GW data when both frequency (f) and its time derivative (ḟ ) are measured. The limiting frequency allowing the chirp mass measurement for a typical BNS is of ∼1.75 mHz (cf. Fig. 2). At lower frequencies, BNS will be seen as monochromatic over the whole duration of the mission, meaning that their chirp mass will be degenerate with the distance (cf. equation 7).
In the case of a circular binary, the measurement of the chirp depends only on f andḟ . The fractional error on the chirp mass can be estimated as where with ρ being the signal-to-noise ratio (SNR) for the mission time T = 4 yr (e.g. Lau et al. 2020). Averaging over sky location, polarization, and inclination, one can write down the SNR as (e.g. Robson, Cornish & Liu 2019) where A is the amplitude of the signal S n (f) is the power spectral density of the detector noise in the low-frequency limit that also accounts for unresolved Galactic background, d is the luminosity distance to the binary, and R(f) is a transfer function encoding finite-arm length effects at high frequencies that computed numerically in Korol et al. (2020). Fig. 1 shows the SNR (cf. equation 6) for a circular BNS placed at the distance of 10 kpc and observed over the nominal 4 yr of the LISA mission. For comparison, we show the subsample of shortest orbital period Galactic BNSs detected through the radio emission (orange stars) with measured masses (Kramer et al. 2006;Ferdman et al. 2014;van Leeuwen et al. 2015;Cameron et al. 2018;Stovall et al. 2018;Farrow et al. 2019). By plugging in their true distances in equation (6), we find that all of the known BNS lay under the LISA detection threshold of ρ = 7. The chirp masses of extragalactic BNS mergers detected through GW emission by LVC are delimited by grey horizontal bands. In Fig. 2, we zoom-in on frequencies >1 mHz and show the expected fractional error on the chirp mass.
The eccentricity is another binary parameter that -when sufficiently high -can impact BNS detectability and, consequently, the chirp mass measurement. In isolated binary evolution, BNSs are expected to form with non-zero eccentricity due to the supernova kicks associated with the formation of the last-born neutron star or Blaauw kicks produced by symmetric mass-loss accompanying supernovae (Blaauw 1961;Tauris et al. 2017). The GW radiation quickly circularizes BNS orbits so that they become almost circular by the time they move to the LISA band (e.g. when the Hulse-Taylor pulsar will evolve to 2 mHz, its eccentricity will decrease from 0.61 to 0.03). However, if the binary is born right in or at the edge of the LISA's frequency window, it will retain the original eccentricity. To assess how the eccentricity influences the chirp mass measurement and, consequently, LISA's ability to constrain the mixing fraction of a bimodal chirp mass distribution, we will explore two limiting cases: the case in which all binaries are circular and the case in which all binaries are eccentric. We anticipate that  Expected fractional error on the chirp mass as a function of GW frequency (x-axis) and chirp mass (y-axis) for a circular binary. Here, we fixed the distance to the binary to d = 10 kpc and the observation time with LISA to T = 4 yr. Overlaid are lines of constant fractional measurement error on the chirp mass; their slope indicates that at a given GW frequency, higher chirp masses lead to a smaller error on the measurement. assuming all binaries to be eccentric does not significantly change our results on the mixing fraction. Thus, we defer the description of how the chirp mass measurement changes in the eccentric case to Appendix A.

Mock Galactic BNS population
To assemble a mock Galactic population, we assume that the BNSs' chirp mass distribution is described by a mixture of two Gaussian distributions: one centred on the value characteristic of the known Galactic population that has been detected at radio wavelengths and another centred on the chirp mass of GW190425. We therefore write where w 1 and w 2 are the relative weights of the two Gaussian distributions normalized such that w 1 + w 2 = 1. We obtain μ 1 = 1.17 M with a standard deviation of σ 1 = 0.04 by combining individual chirp mass measurements of known Galactic BNS reported in Farrow et al. (2019). We set μ 2 = 1.44 M and σ 2 = 0.02 according to the posterior distribution reported in Abbott et al. (2020). We show our model chirp mass distribution in Fig. 3 with the blue solid line. We note that our choice is supported by studies analysing available GW and/or radio observations of BNSs that found evidence for a broad secondary peak at high masses in the birth mass distributions of second-born neutron stars (Farrow et al. 2019;Galaudage et al. 2020).
In addition, population studies of binary white dwarf detectable with LISA also show bimodality in the chirp mass distribution (e.g. Korol et al. 2017).
To model the frequency distribution, we assume that the Galactic BNSs population is stationary on the time-scale of interest. Consequently, their distribution in frequency is given by where R MW is the merger rate of BNSs in the Milky Way. Note, however, that this assumption does not account for possible significant recent star formation episodes that could add BNS systems directly in the LISA band. We can then compute the total number of BNS at frequencies above f by integrating equation (9) N (> f ) = 5c 5 R MW 256π 8/3 (GM c ) 5/3 f 8/3 where we use R MW = 140 Myr −1 (Seto 2019). We note that the Galactic BNS merger rate is still uncertain: based on extragalactic BNS merger rates reported by LVC in the first two observing runs, Andrews et al. (2020)  We note that the updated merger rate of 320 +490 −240 Gpc −3 yr −1 based on the second LVC GW Transient Catalog (The LIGO Scientific Collaboration 2020) corresponds to 32 +49 −24 Myr −1 (assuming the number density of the Milky Way-like galaxies of 0.01 Mpc −3 ), and thus agrees better with the rates estimated from the population of Galactic BNSs. Moreover, in the next decades, the differences in estimated merger rates should further decrease as more observations will become available from both radio and GW observatories.
Finally, we assume that BNS is distributed in the Galactic disc with an exponential radial stellar profile with an isothermal vertical distribution where 0 kpc ≤ R ≤ 20 kpc is the cylindrical radius measured from the Galactic Centre, z is the height above the Galactic plane, R d = 2.5 kpc is the characteristic scale radius, and z d = 0.4 kpc is the vertical scale height of the observed Galactic BNSs (Pol et al. 2019). Finally, when converting BNS positions (R, z) into heliocentric distances (d), we assume the position of the Sun to be at (8.1,0) according to Gravity Collaboration (2019).

Bayesian inference
Given N measurements of the BNS chirp masses with associated errors, we now would like to reconstruct the shape of the underlying true chirp mass distribution. Here, we have chosen to model the chirp mass of Galactic BNSs as a mixture of two Gaussian distributions (cf. equation 8), thus the distribution is fully described by six parameters λ ∈ {w 1 , w 2 , μ 1 , μ 2 , σ 1 , σ 2 }. We follow the Bayesian approach as outlined in Mandel, Farr & Gair (2019) ignoring selection effects because we are only interested in a subsample of detections with f > 2 mHz -that allows the measurement of the chirp massacross which we find the detection probability with LISA to be ∼1.
To simulate the instrumental noise, we displace each chirp mass from its true value by resampling it from a Gaussian centred on the true value and standard deviation determined by the LISA's measurement error σ M (cf. equation 3). We will denote this displaced chirp mass withM.
Using Bayes' theorem, we can write the posterior as where π (λ) are priors on {w 1 , w 2 , μ 1 , μ 2 , σ 1 , σ 2 }, P (M i |M) is the probability of observing the event i given the assumed underlying distribution (likelihood), and P pop is the probability of individual chirp mass given the underlying population distribution (equation 8).

R E S U LT S
We start by assuming the fraction of BNSs that resemble GW190425 to be w 2 = 0.2 and we set the Galactic merger rate to be 140 Myr −1 , which corresponds to 33 BNSs with frequencies higher than 2 mHz in the Galaxy (cf. equation 10). We recover the parameters of the chirp mass distribution in the framework of probabilistic programming package PYMC3 (Salvatier, Wiecki & Fonnesbeck 2016) and sample the posterior (equation 12) using the No U-turn Sampler. Fig. 4 shows two examples of posterior distributions for w 2 , μ 1 , μ 2 , σ 1 , and σ 2 . We have chosen to display only one of the two weights w 2 , since the other can be recovered as 1 − w 2 . We present our results for two assumptions on the BNS merger rate: R MW = 42 Myr −1 in orange and R MW = 212 Myr −1 in blue colours. The true values are marked by black dashed lines. It is immediately evident that the higher merger rate leads to more constrained results because the number of detectable BNSs by LISA is higher. Moreover, in the case of the lower merger rate, the small number statistics leads to bias in some of the parameters, although the true value is recovered to within the 1σ of the derived posteriors. Importantly, Fig. 4 shows that we can constrain the relative weights of the two BNS subpopulations w 1 and w 2 . We now focus on how well our inference machinery can recover w 2 (or equivalently w 1 ). We fix the Galactic merger rate and vary w 2 . We find that for R MW = 140 Myr −1 , the weight of GW190425like subpopulation is recovered with 1σ error of 0.1-0.04, where the largest error is obtained for w 2 = 0.1, while the smallest for w 2 = 0.9. This trend can be explained by remembering that the error on the chirp mass is mainly dominated by σḟ ∝ M −5/3 (cf. equation 5 for the circular case; for the eccentric case this is true only up to ∼3 mHz, where σḟ and σ e start to be comparable, see Appendix A). Thus, the population with w 2 = 0.9 has overall better measured chirp masses than a population with w 2 = 0.1. Next, we fix w 2 and vary the merger rate in the range between 21 and 212 Myr −1 . We find that the error on the recovered fraction decreases with increasing merger rate as more LISA observations become available. Specifically, when setting the 'true' value w 2 = 0.2, we recover 0.34 +0. The summary of the recovered w 2 as a function of the 'true' (input) w 2 is represented in Fig. 5. It shows that chirp mass measurement errors allow the recovery of the two subpopulations' weights with no bias regardless of its 'true' value. In colour, we show our ability to recover w 2 for different Galactic merger rates of 42, 140, and 212 Myr −1 , which corresponds to the total number of observed binaries at f > 2 mHz of 10, 33, and 50, respectively (cf. equation 10).
Finally, for comparison, we perform the same set of simulations as presented above for the limiting case in which all BNSs in the LISA band are eccentric. For example, using different natal kick prescriptions, Lau et al. (2020) report median eccentricity of the population ranging from 0.36 to 0.071, with a median of 0.1 for their fiducial model. For simplicity, as an example here, we set the eccentricity of all binaries to 0.3. We do not find any significant deviations in recovering w 2 from the results of the circular case presented above. We can attribute this to the fact that the chirp mass error distribution does not differ significantly between the circular and eccentric cases (see also Appendix A).

D I S C U S S I O N A N D C O N C L U S I O N S
In this paper, we have investigated whether future observations of Galactic BNS with LISA and specifically LISA's chirp mass measurement can elucidate the nature of heavy BNS like GW190425. First, we showed that if GW190425-like binaries populate frequencies >2 mHz, they can be detected by LISA within the volume of our Galaxy. Moreover, their chirp masses can be measured to better than 10 per cent. Then, we constructed a toy model for Galactic BNSs consisting of two distinct subpopulations: one that follows the chirp mass distribution of known BNSs detected through radio observations and another that resembles GW190425. Since the relative fraction of the two subpopulations depends on the BNS formation model and the treatment of the fast-merging channel in modelling the total population in the Galaxy (e.g. Vigna-Gómez et al. 2018;Mandel et al. 2021), here we considered the weights of the two subpopulations as free parameters. As the Galactic BNS merger rate is still uncertain, we repeat our study for a range of merger rates between 24 and 212 Myr −1 quoted in the literature (Abbott et al. 2017;Vigna-Gómez et al. 2018;Pol et al. 2019;Andrews et al. 2020). We demonstrated that if GW190425-like binaries constitute a fraction larger than 0.1 of the total population, LISA should be able to recover the fraction with better than ∼15 per cent accuracy assuming the merger rate of R MW = 42 Myr −1 (corresponding to 10 detected binaries with f > 2 mHz); the accuracy increases to ∼5 per cent for R MW = 212 Myr −1 (corresponding to 50 detected binaries with f > 2 mHz). We note that with more upcoming radio and GW observations, the BNS merger rate is expected to be better constrained, such that we will have a more precise estimate on the expected number of BNS sources detectable by LISA. The results of this work can then be used to access what fractions of more massive GW190425-like binaries can be constrained by LISA. Finally, our results show that even if all Galactic BNSs come with large eccentricities, the errors on the measured chirp masses stay within the limits where our recovered fractions, assuming circularized orbits, remain applicable (cf. Fig. 6).
Further considerations on BNS formation could be deduced from the comparison of the BNS rate determined from radio observations of pulsars at sub-mHz frequencies (stars in Fig. 1) with that inferred from LISA at mHz frequencies. If all BNSs form at low frequencies (merger times of >100 Myr) and evolve to mHz frequencies via GW emission (cf. equation 10), the two rates should be consistent. However, if the BNS rate at mHz frequencies is larger than that predicted by equation (10), the excess could be attributed either to a subpopulation of radio-quiet BNSs (e.g. Safarzadeh et al. 2020) or to the fast-merging channel.
Population synthesis studies suggest a broad distribution of delay times between BNS formation and merger even under the assumption that the case-BB mass transfer is dynamically stable (e.g. Tauris et al. 2013;Tauris, Langer & Podsiadlowski 2015;Vigna-Gómez et al. 2018). However, if the case-BB mass transfer is occasionally dynamically unstable (Dewi & Pols 2003;Ivanova et al. 2003), it is possible that fast-merging BNSs form within the LISA band at higher frequencies than the rest of the Galactic population. Therefore, if their merger time-scales are sufficiently long to be observed with LISA, it is possible that fast-merging BNSs exhibit an excess in the frequency distribution. LISA's frequency measurement is expected to be very accurate (ideally σ f /f ∝ 1/T ∼ 10 −8 for T = 1 yr); thus, identifying a pile-up of BNSs at high frequencies should be a straightforward task. Importantly, the observed frequency distribution will facilitate constraints on the BNS delay time distribution. Moreover, highfrequency binaries would also have better chirp mass measurements (and other parameters including the eccentricity) that could only improve the results presented here.
To confidently constrain the fast-merging hypothesis, eccentricity measurements are required (Andrews & Mandel 2019;Andrews et al. 2020;Lau et al. 2020). If high-mass GW190425-like binaries are also highly eccentric (e ∼ 1), this would imply that they are freshly formed (either in isolation or in a dense cluster), and thus come from the fastmerging channel. Fig. 6 shows that constraining both the chirp mass and the eccentricity to, for example, 10 per cent is only possible at mHz frequencies. Andrews et al. (2020) showed that such fractional errors on the eccentricity should be sufficient to distinguish between the different BNS formation channels. Still, if the eccentricity of GW190425-like binaries is low (e ∼ 0), one cannot rule out the possibility that they are a population of radio-quiet BNSs.
We remark that our error estimates are based on the Fisher information matrix approximation and are valid for high SNRs (e.g. Cutler 1998). Thus, we expect that in some cases the reported errors may be underestimated. A full Bayesian parameter estimation would be required to derive more realistic uncertainties (e.g. Buscicchio et al. 2019;Roebber et al. 2020).
In this work, we focused exclusively on the LISA mission, but our analysis would hold also for the TianQin (Luo et al. 2016;Mei et al. 2020) and Taiji (Ruan et al. 2018) space-based GW observatories designed to operate in a similar frequency range as LISA. Huang et al. (2020) showed that LISA and TianQin can simultaneously detect Galactic stellar binaries above a few mHz and that combined observations from the two missions will improve their parameter estimation including orbital period, inclination, and sky localization.
While radio observatories are sensitive to the intermediate stages of BNS evolution and ground-based GW detectors can see only the very last seconds of a BNS's life and the final merger, LISA will be crucial for bridging the gap between these two regimes (Galaudage et al. 2020). Here, we argued that the BNS chirp mass distribution measured by LISA could be a useful tool for unveiling massive GW190425-like BNSs in the Milky Way and constraining their origin. Other studies established the importance of the LISA's eccentricity measurement for understanding the BNS formation pathways (Andrews et al. 2020;Lau et al. 2020). In addition, LISA will also provide sky positions for Galactic BNS enabling radio follow-up and the discovery of radio-faint pulsars with all-sky surveys such as the Square Kilometre Array (Kyutoku, Nishino & Seto 2019). Finally, the contribution to chirp mass error due to eccentricity σ F(e) can be calculated directly for known σ e using equation (A2).
In the top panel of Fig. A1, we show how the fractional error on the chirp mass changes as a function of the eccentricity for a binary at the distance of d = 10 kpc and emitting at 2 mHz. In the bottom panel of Fig. A1, we fix binary's eccentricity to e = 0.3 and split fractional error of the chirp mass (black) into the contributions due toḟ (blue) and e (orange); we do not represent the contribution due to the uncertainty as it remains negligible at all frequencies. For f 3 mHz, the error on the chirp mass is dominated by the error onḟ , which scales with ρ and thus σ M /M is smaller compared to in the circular case (see upper panel). At 3-4 mHz, the contribution of the eccentricity limits the chirp mass measurement, so at higher frequencies the chirp mass is better for the circular binaries.