{'fraction_non_alphanumeric': 0.03542034993598732, 'fraction_numerical': 0.02341035176492105, 'mean_word_length': 4.312176165803109, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constraining the Milky Way potential with Double White Dwarfs

The upcoming LISA mission is the only experiment that will allow us to study the Milky Way's structure using gravitational wave signals from Galactic double white dwarfs (DWDs). The total number of expected detections exceeds $10^5$. Furthermore, up to a hundred DWDs can be simultaneously detected in both gravitational and optical radiation (e.g. with Gaia and LSST as eclipsing), making DWDs ideal sources for performing a multi-messenger tomography of the Galaxy. We show that LISA will detect DWDs everywhere, mapping also the opposite side of the Galaxy. This complete coverage will: (1) provide precise and unbiased constraints on the scale radii of the Milky Way's bulge and disc, and (2) allow us to compute the rotation curve and derive competitive estimates for the bulge and disc masses, when combining gravitational wave and optical observations.


Introduction
So far we have only explored the Milky Way through electromagnetic (EM) radiation, but the upcoming LISA mission will allow us to probe our Galaxy with gravitational wave (GWs) radiation. Moreover, by exploiting GWs we can circumvent some difficulties inherent to EM observations, such as the dust extinction at optical wavelengths, luminosity bias and various difficulties related to distance determination. LISA is a space-based GW experiment approved by ESA in 2017 and scheduled for launch in early 2030 (Amaro-Seoane et al., 2017). Its design consists of a set of three identical spacecrafts in an equilateral triangle formation, with a center of mass that will follow a circular heliocentric trajectory lagging behind the Earth. The spacecraft separation is set to 2.5 × 10 6 km, making LISA sensitive to GW frequencies between 0.1 and 100 mHz. Our Galaxy hosts a variety of stellar mass mHz GW sources like detached DWD binaries, AM CVn systems, ultra-compact X-ray binaries and binary black holes. A number of these binaries have been identified though EM radiation (e.g., Kupfer et al., 2018). In particular, detached DWDs are foreseen to be the most abundant LISA sources, with > 10 5 individually resolvable in frequency (e.g., Nelemans et al., 2004;Ruiter et al., 2010;Marsh, 2011;Korol et al., 2017). Indeed, these binaries are predicted to be so common in the Milky Way that overlapping signals of those below the LISA sensitivity threshold will form a background noise for the mission (e.g., Robson & Cornish, 2017). Both resolved DWDs and the DWD background bear the imprint of the overall Galactic stellar population, thus both can be used to study the Milky Way's structure and star formation history. By fitting density distributions of resolved LISA detections we can precisely estimate the scale radii of the disc and the bulge, while by fitting the DWD background we can recover the disc scale height (Benacquista & Holley-Bockelmann, 2006;Adams et al., 2012;Korol et al., 2018). Additional information (such as the motion of DWDs) required to constrain the Milky Way potential can be recovered from optical observations. By fitting the Galactic rotation curve constructed using distances inferred from gravitational waves and proper motions from optical observations one can obtain competitive estimates of the bulge and the disc masses .

Synthetic Milky Way
To obtain a mock DWD population we exploit the population synthesis code SeBa (Portegies Zwart & Verbunt, 1996;Toonen et al., 2012). The initial stellar population is obtained assuming the Kroupa initial mass function, a flat binary mass ratio distribution, a log-flat distribution for the binary semi-major axis, a thermal distribution for the orbit eccentricity, an isotropic distribution for binary inclination angles, and a constant binary fraction of 50% (Kroupa et al., 1993;Duchêne & Kraus, 2013;Raghavan et al., 2010). We use the γα prescription for the common envelope phase (Nelemans et al., 2001). The sensitivity of our model to these assumptions is discussed in Korol et al. (2017) and Toonen et al. (2017). We adopt a potential composed of a bulge, a stellar disc and a dark matter (DM) halo. The density distribution of the disc component 2 Constraining the Milky Way potential with Double White Dwarfs for our model can be analytically expressed as (1) where 0 ≤ R ≤ 19 kpc is the cylindrical radius measured from the Galactic center, R d = 2.5 kpc is the characteristic scale radius, Z is the height above the Galactic plane and Z d = 300 pc is the characteristic scale height of the disc. To account for the star formation history we use the plane-projected star formation rate from Boissier & Prantzos (1999), ρ BP , assuming the age of the Galaxy to be 13.5 Gyr. The total mass of the disc in our model is 5 × 10 10 M and the distance of the Sun from the Galactic center is R = 8.5 kpc. The bulge component is modeled by doubling the star formation rate in the inner 3 kpc of the Galaxy and distributing DWDs spherically: where r is the spherical distance from the Galactic center, M b = 2.6 × 10 10 M is the bulge total mass, and r b = 0.5 kpc is the characteristic radius of the bulge. Finally, to describe the DM halo we use the Navarro et al. (1996) density profile: where r s = 20 kpc is the scale length of the halo and ρ h = 0.5 × 10 7 M kpc −3 is the halo scale density. The total mass of the halo enclosed in a 100 kpc radius is equal to 4.8×10 11 M . We summarize the values of the parameters adopted for our Milky Way fiducial model in Table 1.

Milky Way map in gravitational waves
To model GW signals of DWDs we employ the Mock LISA Data Challenge (MLDC) pipeline designed for the simulation and analysis of a large number of quasimonochromatic GW sources (e.g., Littenberg, 2011). The MLDC pipeline describes GW signals in terms of nine GW parameters: amplitude A, frequency f = 2/P with P being the binary orbital period,ḟ ,f , sky location in the LISA reference frame (θ, and φ), orbital inclination ι, GW polarization ψ and the initial phase φ 0 . For a given synthetic instrument noise curve, observation time and detection threshold, the pipeline provides a catalog of sources with signal-to-noise ratio (SNR) above the detection threshold and the uncertainties of the nine GW parameters. We adopt a mission duration of 4 yr and the noise curve from Amaro-Seoane et al. (2017), which corresponds to the mission design approved by ESA. Finally, we set the SNR detection threshold to 7. We find 2.6×10 4 DWDs in our mock Milky Way above the threshold. In particular, we estimate their distances and sky localization with their respective uncertainties. We compute the distances from GW observables A, f andḟ as (Maggiore, 2008) and the respective errors as where σ A /A, σ f /f and σḟ /ḟ are outputs of the MLDC pipeline. We verify that the terms containing correlation coefficients are at most of the order of 1%, and we thus neglect them in eqn. (5). The coordinate transformation between the LISA reference frame (θ, φ) and the Galactic coordinate frame (l, b) is outlined in Korol et al. (2018). The distribution of LISA detections in Galactic coordinates is shown in Fig. 1. The map reveals a spherical agglomeration of stars in the center, which represents the bulge, and an extended flat distribution of stars, which represents the stellar disc. We find that LISA will be able to detect these binaries to large distances, comparable with the extension of the disc, also mapping of the opposite side of the Milky Way (see Korol et al., 2018, figure 3). Out of 2.6 × 10 4 LISA detections in our simulation, almost 8×10 3 DWDs (30%) have relative distance errors < 30%. In particular, a subsample of ∼ 100 DWDs (0.4% of the catalog) has relative errors on the distance < 1%. These are high frequency (> 3 mHz), high SNR (> 100) sources spread around the Galaxy. This remarkable precision is due to the fact that GW amplitudes decrease more slowly with distance (∝ 1/d) than EM observations (∝ 1/d 2 ). This property of GWs is at the heart of the unique ability of the LISA mission to study the Milky Way's structure. Finally, we find that on average LISA will locate DWDs within 1 − 10 deg 2 on the sky. Figure 1 suggests that LISA has the potential to reconstruct the density profiles of both the disc and bulge components of our mock Galaxy. To test this we select binaries obtained from the LISA catalog with distances determined to better than 30%. First, we derive the radial density profile. We compute 10 5 realizations of the 3D binary positions in the Galaxy by randomly drawing l, b and d from Gaussian distributions centered on their true values and with standard deviations provided by the MLDC pipeline. For each realization we compute the cylindrical Galactocentric distance, R, and we select those with 2 ≤ R ≤ 12 kpc. We bin DWDs in cylindrical shells of width dR = 125 pc, and divide the bin counts by the shell volume. We compute the error on the number density in each bin as the standard deviation over different realizations. The obtained disc radial density profile is shown in the upper panel of Fig. 2. Then, we fit the density profile with the analytic expression in eqn.

Density distributions
(1), obtaining the value for the disc scale radius of R d = 2.54±0.08 kpc, in agreement with the fiducial value. The blue solid curve shows the best fit, and the colored area shows the 3σ interval.
Next, we study the vertical distribution of DWDs in the disc. Again, we select only binaries with 2 ≤ R ≤ 12 kpc, and we bin them in the R − Z plane. In each radial bin, we model the number density with a sech 2 (Z/Z d ) function and fit Z d to test whether the scale height is constant with R. We find a constant behavior and therefore we decide to increase the statistics by computing the average value of Z d and its error by using a stacked radial profile. In this way, we find Z d = 0.31 ± 0.04 kpc, which is consistent with the fiducial value of 0.3 kpc.
Finally, to estimate the scale radius of the bulge we select DWDs in the inner 1 kpc, and we compute 10 5 realization of the binary positions in the Galaxy, as we did for the disc radial profile. For each realization we estimate the number density profile by counting DWDs in spherical shells with radius r and dr = 20 pc, dividing this number by the shell volume. The result is given by the magenta triangles in the bottom panel of Fig. 2. To fit the scale radius of the bulge, we use eq. (2) as the model distribution, and we obtain r b = 0.51 ± 0.013 kpc.

Combining gravitational wave and optical observations
The mass of the Milky Way's disc and bulge cannot be derived from GW data alone. Optical observations of the motion of DWDs in the sky are required to constrain them. The sky localization from GW data is typically poor, which makes it difficult to identify counterparts to GW sources in EM databases. In practice, to assemble a sample of optical counterparts, one possibility is to search in optical surveys for periodically variable sources with a frequency and within an area on the sky matching those provided by LISA. To assess whether this is possible we focus on edge-on binaries, which allow for better distance estimation with GWs and are also easy to identify by optical telescopes as eclipsing sources (Shah et al., 2012). In our previous study we found that Gaia (Gaia Collaboration et al., 2016) and the Large Synoptic Survey Telescope (LSST, LSST Science Collaboration et al., 2009) can identify 75 of DWDs detected by LISA as eclipsing (of which 23 are identified by both) (Korol et al., 2017). We use this sample of binaries also in this work. To understand how well it will possible to constraint the motion of the DWD optical counterparts we simulate their parallaxes and proper motions µ as measured by Gaia and LSST. We draw the parallaxes from a Gaussian distribution centered on 1/d and with standard deviation σ . We 4 Constraining the Milky Way potential with Double White Dwarfs  compute proper motions as µ = V t /d, where V t is the tangential velocity, which for a fixed Galactic potential can be derived from geometrical consideration as outlined in Korol et al. (2018). The end-of-mission parallax and proper motion errors (σ and σ µ ) for Gaia can be estimated with the pyGaia 1 python toolkit, and for LSST they are tabulated in LSST Science Collaboration et al. (2009). For those DWDs that are identified by both optical instruments we select the measurement with the smallest uncertainty. Next, we convert parallaxes into distances following the Bayesian inference approach indicated in Bailer-Jones (2015); Luri et al. (2018), adopting an exponentially decreasing space density function with scale length of 400 pc as the prior.
Our choice for the scale length is explained in (Kupfer et al., 2018). We find that the expected relative error in parallax for binaries at d < 1 kpc is < 20%, however this is only 30% of the optical counterparts. This gives us the motivation to combine distances derived from parallaxes with the LISA measurements. Again, this can be done by using Bayes' theorem and is described in Korol et al. (2018), Sect. 3.3. The resulting distances (hereafter d obs ) represent the harmonic means of the two measurements (GW and optical) and the respective errors are equal to twice the harmonic mean of the individual errors. Essentially, in this way we select the best between GW and optical measurements and we reduce the uncertainty on the distance compared to just selecting the more precise of the EM or GW measurements individually.

Rotation curve
To simulate the motion of DWD in the Galaxy we compute the observed rotation speed (in km/s) as  motions, we assign an observed proper motion µ obs to a source by sampling from a Gaussian centered on µ and with an error σ µ . For each DWD, we calculate V obs (R) for 10 5 independent realizations of µ obs and d obs , and we assign an observed velocity and a measurement error equal, respectively, to the mean and the standard deviation of the resulting distribution of V obs (R). The result is represented in Fig. 3. We fit the obtained data with the total rotation curve corresponding to our Milky Way model, which we compute as V 2 (R) = RdΦ/dR with Φ being the total potential. We perform the fit with an MCMC code fixing r b , R d and Z d to the values obtained in Sect.3.1, and leaving M d and M b , ρ h , and r h free, but confined by flat uninformative priors in the respective ranges: (1 − 10) × 10 10 M , (0.1−10)×10 7 M /kpc 3 and 10−30 kpc. The final posterior probability distribution for the free parameters is represented in Fig. 4. It reveals that we can recover the mass of the disc and bulge components, but not that of the DM halo. This is because there is no optical data at R > 11 kpc, where the halo component dominates the dynamics in our Milky Way model (Fig.3). We estimate the mass of the disc to be M d = 5.3 +1.29 −1.71 ×10 10 M and the mass of the bulge to be M b = 2.49 +0.44 −0.42 × 10 10 M , in good agreement with our fiducial values. Remarkably, our constraints on the bulge mass are extremely competitive with those derived from EM tracers (see e.g. Bland-Hawthorn & Gerhard, 2016, for a review). The larger errors on the disc mass stem from our choice to leave the halo parameters unconstrained. Thus, an improvement of this analysis should involve the adoption of stringent constraints on the DM halo parameters. We summarize all the recovered parameters in Tab. 1.

Conclusions
In this study we introduce the idea of the multimessenger study of the Milky Way using Galactic DWD binaries, and we investigate the prospects for tracing the baryonic mass of the Milky Way using GWs in combination with their optical counterparts. The advantages over traditional tracers include the possibility of looking through the bulge, and beyond, thus allowing one to map both sides of the Galaxy using the same tracer. We show that the unique property of DWDs allows one to recover the scale radii of the baryonic components with a percent precision. In synergy with optical data, GW measurements will provide competitive mass estimates for the bulge and stellar disc. However, our choice to use GW sources and their EM counterparts limits our ability to constrain the DM halo component of the Galaxy. This highlights the importance of a more precise knowledge of the DM halo to improve baryonic mass measurements.