The Milky Way's bar structural properties from gravitational waves

The Laser Interferometer Space Antenna (LISA) will enable Galactic gravitational wave (GW) astronomy by individually resolving $>10^4$ signals from double white dwarf (DWD) binaries throughout the Milky Way. In this work we assess for the first time the potential of LISA data to map the Galactic stellar bar and spiral arms, since GWs are unaffected by stellar crowding and dust extinction unlike optical observations of the bulge region. To achieve this goal we combine a realistic population of Galactic DWDs with a high-resolution N-Body simulation a galaxy in good agreement with the Milky Way. We then model GW signals from our synthetic DWD population and reconstruct the structure of the simulated Galaxy from mock LISA observations. Our results show that while the low signal contrast between the background disc and the spiral arms hampers our ability to characterise the spiral structure, the stellar bar will instead clearly appear in the GW map of the bulge. The bar length and bar width derived from these synthetic observations are underestimated, respectively within $1\sigma$ and at a level greater than $2\sigma$, but the resulting axis ratio agrees to well within $1\sigma$, while the viewing angle is recovered to within one degree. These are competitive constraints compared to those from electromagnetic tracers, and they are obtained with a completely independent method. We therefore foresee that the synergistic use of GWs and electromagnetic tracers will be a powerful strategy to map the bar and the bulge of the Milky Way.


INTRODUCTION
It is established that the Milky Way has a central stellar bar and a spiral structure that propagates through its stellar and gaseous disc. Observations of molecular masers associated with very young high-mass stars strongly suggest that the Milky Way is a four-arm spiral (Reid et al. 2019). However the nature and the structure of these non-axisymmetric features remain uncertain. Indeed, the amplitude, length and pattern speed of the stellar bar are debated, and constant effort to determine their values is motivated by their importance for a broad range of Galactic studies. For example, these bar features are key to understanding the properties of the disc outside the stellar bar (Minchev & Famaey 2010), the kinematics in the solar neighbourhood (e.g., Dehnen 2000;Pérez-Villegas et al. 2017;D'Onghia & Aguerri 2019;Monari et al. 2019), and the observed non-circular gas flow (e.g., Bissantz et al. 2003).
Thanks to the Gaia mission -(Gaia Collaboration et al. 2016Collaboration et al. , 2018 that is releasing proper motions and distances for almost two billion stars -progress has been made in understanding the central region of the Galaxy, despite the strong dust extinction (e.g., Anders et al. 2019;Bovy et al. 2019). The current picture of the central Contact e-mail: wilhelm@strw.leidenuniv.nl Galactic region is that the stellar bar extends as far as 4.5 − 5 kpc from the Galactic centre and the pseudo-bulge is likely peanutshaped (Benjamin et al. 2005;Nataf et al. 2010;Zoccali & Valenti 2016;Wegg et al. 2015;Valenti et al. 2016).
Of great importance for the dynamics of the bar and its surrounding disc is the pattern speed Ω P of the bar itself (or equivalently its co-rotation radius), because it defines the locations of resonances associated to the bar (e.g., D'Onghia & Aguerri 2019, table 1). Past measurements indicated a rather fast rotation rate (pattern speed) for the Galactic bar of the order of Ω P = 55 km s −1 kpc −1 (Englmaier & Gerhard 1999;Fux 1999;Debattista et al. 2002;Bissantz et al. 2003) with a bar length of ∼ 3 kpc that suggested the Sun being close to the outer Lindblad resonance of the bar.
However, recent measurements of both the three-dimensional density of red clump giants (Wegg et al. 2015) and the gas kinematics in the inner Galaxy (Sormani & Magorrian 2015) point out that the bar pattern speed might be significantly slower. In particular, Portail et al. (2017), by modelling the kinematics and photometry of stars in the bar, inferred a pattern speed of 39 km s −1 kpc −1 . Recently, Sanders et al. (2019) derived a pattern speed of Ω P ≈ 41 km s −1 kpc −1 , applying a modified Tremaine & Weinberg (1984) method on a proper motion data set assembled from both multi-epoch VVV survey and Gaia data release 2 (DR2). Another study derived line-of-sight integrated and distanceresolved maps of average proper motions and velocity dispersions from the VVV Infrared Astrometric Catalogue, combined with data from Gaia DR2, and found agreement with a bar pattern speed of 39 km s −1 kpc −1 (Wegg et al. 2015). A combination of Gaia DR2 and APOGEE data has been also used to infer a similar bar pattern speed (Bovy et al. 2019), although measurements may be affected by systematic uncertainties resulting from photometric distances and the small number of Gaia stars. The Bayesian code S H , combining Gaia parallax information with APOGEE spectroscopic information in addition to other photometric bands, obtained distances in the region dominated by the bar and bulge (Anders et al. 2019). Yet, distances are inferred with uncertainties of the order of 1 kpc (∼10 per cent precision).
Another unresolved issue concerns the bar 'viewing angle', i.e. the orientation of the Sun with respect to the long-axis of the bar. Most authors have found a viewing angle in the range between 10 -45 • (Simion et al. 2017, figure 17).
Recently, it has been shown that it is possible to re-construct a full 3D picture of the Milky Way by exploiting gravitational wave (GW) radiation from Galactic ultra-compact detached double white dwarf (DWD) binaries (Korol et al. 2019). Galactic DWDs can be detected in the milli-Hz GW band with the future Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017). Population synthesis studies forecast > 10 4 DWDs to be individually resolvable by LISA (e.g. Korol et al. 2017;Lamberts et al. 2019;Breivik et al. 2019a). Such a large number of detections distributed across the Galaxy will allow us to map the Milky Way in GWs and precisely measure its structural parameters like the scale radii of the bulge and the disc (Adams et al. 2012;Korol et al. 2019), while combined GW and optical observations of DWDs can be used to derive the mass of the bulge and the disc component of the Galaxy (Korol et al. 2019). At frequencies <3 mHz the LISA band starts to be overpopulated by DWDs giving rise to an unresolved confusion background (e.g., Robson & Cornish 2017). Although affecting detectability of extragalactic LISA sources, this background encodes the properties of the overall stellar population in the Galaxy and can also be used to recover Milky Way's parameters such as the disc scale height (Benacquista & Holley-Bockelmann 2006;Breivik et al. 2019b).
In this work we explore for the first time the potential of GWs to characterise the Milky Way's bar and spirals structural properties. To achieve this we combine a DWD binary population synthesis model with GALAKOS, a high resolution Milky Way-like galaxy simulation (D'Onghia & Aguerri 2019). We consider the population of detached DWDs only, because distance determination for these binaries is not affected by systematic uncertainties due to astrophysical processes such as mass transfer. In Section 2 we explain how we combine these two tools to obtain a realistic catalogue of Galactic DWDs. In Section 3 we describe how the catalogue is processed to obtain mock LISA observations. Next, we outline the method adopted to infer the structure of the Galaxy from LISA observations (Section 4). In Section 5 we present our results. Finally, we discuss how our results compare to electromagnetic observations and draw conclusions in Section 6.

METHODS
To model the population of DWDs in the Milky Way we follow the same method as in Korol et al. (2019), but instead of employing an analytic gravitational potential we use high resolution numerical simulations of D'Onghia & Aguerri (2019). This Galaxy model  (Boissier & Prantzos 1999). The present day value is 1.87 M yr −1 .
includes a self-consistent spiral pattern and bar whose properties are in agreement with current observations. It allows us to test how the Galactic structure will be observed with LISA.

Binary population synthesis model
We employ the DWD population model from Toonen et al. (2012), obtained using binary population synthesis code S B , developed by Portegies Zwart & Verbunt (1996), later adapted for DWDs by Nelemans et al. (2001) and Toonen et al. (2012). The progenitor population is constructed by randomly sampling binary properties with a Monte Carlo technique from distributions motivated by currently available observations for intermediate-and low-mass stars. Specifically, the mass of the primary star is drawn from the initial mass function of Kroupa et al. (1993) in the range between 0.95 and 10 M . The mass of the secondary star is derived from a uniform mass ratio distribution between 0 and 1 (Duchêne & Kraus 2013).
To draw the initial binary orbital separations and binary eccentricities we adopt respectively a log-flat and a thermal distributions (Abt 1983;Heggie 1975;Duchêne & Kraus 2013). We set the initial binary fraction to 0.5. First, S B evolves the initial population from Zero-age mainsequence until both stars become white dwarfs. Then, from DWD formation until the present time binaries are further evolved via GW emission that causes binaries' orbits to shrink (e.g., Peters 1964). Finally, we remove binaries from the catalogue if they have begun mass transfer (i.e. when one of the two white dwarfs fills its Roche lobe) or they have already merged within the present time.
To assign a realistic present time age distribution of DWDs we used the star formation rate (SFR) grid from a hydrodynamic simulation of Boissier & Prantzos (1999). Our modelling of star formation implicitly assumes that once the bar is formed, it drives gas toward the central regions of the Galaxy enhancing star formation. This is a reasonable assumption as the majority of the bulge forms as part of the bar, according to the recent estimates of the kinematics of the stars of the bulge and bar ). Our method does not envisage a different star formation history in the Milky Ways's spiral arms. However, as the bar and spiral structure develop the stars orbiting in the disc tend to increase the amplitude of their epicycle causing an increase of the in-plane velocity dispersion and a global heating of the disc (Binney & Tremaine 1987). As a result the spiral arms that propagate through the old stellar disc formed also by the old stellar population such as DWDs will be less visible.
For the distribution of DWDs in our Milky Way model we resampled the DWD population of Korol et al. (2019), without information whether a DWD originated in the bulge or in the disc. Thus, our population has a single star formation history equal to the (weighted) sum of the disc and bulge star formation histories of Boissier & Prantzos (1999), for both its disc and bulge.
One of the most impacting assumption in binary population synthesis is the prescription for the common envelope (CE) evolution. CE is a short phase of the binary evolution (∼ 10 3 yr) in which the more massive star of the pair expands and engulfs its companion (Paczynski 1976;Webbink 1984). During this phase the binary orbital energy and angular momentum can be transferred to the envelope, due to the dynamical friction experienced by the companion star experiences when moving through the envelope. This process continues until the envelope is ejected from the system leaving behind the core of the expanded star and its companion in a tighter orbit. Typically, the CE is implemented in the binary population synthesis either by parametrising the conservation equation for energy (through the α parameter) or the balance equation of the system's angular momentum (through the γ parameter) (see Ivanova et al. 2013, for a review). For this study we adopt the γα DWD evolution model, in which both parametrisations are allowed: the γ-prescription is applied unless the binary contains a compact object or the CE is triggered by a tidal instability, in which case α-prescription is used. Synthetic catalogues of DWDs produced using this evolution model have been carefully calibrated against state-of-the-art observations of DWDs in terms of both mass ratio distribution (Toonen et al. 2012) and number density . Future optical surveys such as the Vera Rubin Observatory (Ivezic et al. 2008) will provided large samples of new DWDs that will help to further constrain DWD evolution models (Korol et al. 2017).

Milky Way model
To simulate the Milky Way we use a snapshot of GALAKOS, a high-resolution N-body simulation of a stellar disc with structural parameters that reproduce the currently observed properties of our Galaxy (D'Onghia & Aguerri 2019). The simulation was carried out with GADGET3, a parallel TreePM-Smoothed particle hydrodynamics (SPH) code developed to compute the evolution of stars and dark matter, treated as collisionless fluids. The phase space is discretised into fluid elements that are computationally realised as particles in the simulation. The total number of N-Body particles employed in the simulation is 90 million. The gravitational softening length adopted is 40 pc for the dark halo, 28 pc for the stellar disk and 80 pc for the bulge. Note that the gas component is not included in the simulation. The Milky Way model consists of three components: a dark matter halo, a rotationally supported stellar disc, and a spherical stellar bulge. These components are modelled as follows: • Dark matter halo is modelled with the Hernquist (1990) density profile where M DM = 1 × 10 12 M is the total dark matter mass and a = 30 kpc is the radial scale length. This choice is motivated by the fact that in its inner part, the shape of the density profile is identical to the Navarro-Frenk-White fitting formula of the mass density distribution of dark matter halos inferred in cosmological simulations (Navarro et al. 1996). The dark matter mass in GALAKOS is sampled with 60 million particles.
• Stellar disc is represented by an exponential radial stellar disc profile with an isothermal vertical distribution where M d = 4.8 × 10 10 M is the disc mass, R d = 2.67 kpc is the disc scale length, and z 0 = 320 pc is the disc scale height. The disc mass is discretised with 24 million particles. The gas component is not included in the model.
• Stellar bulge is also described by the Hernquist model. The total mass of the bulge adopted in the simulation is M b = 8 × 10 9 M with a scale length a b = 320 pc. The number of particles sampled in the bulge is 8.4 million. Note, that the bulge in this simulation does not rotate.
The simulation results in a Milky Way-like galaxy with the total stellar mass of 5.6 × 10 10 M and accounts for a time-varying potential that after 2.5 Gyrs forms a bar with a length of 4.5 kpc, a width of ≈ 2.5 kpc and a pattern speed of 40 km s −1 kpc 1 (see Appendix in D'Onghia & Aguerri (2019) for the details). Figure 2 shows the Milky Way model under different projections. Note that this model includes two prominent spiral arms, whereas the Milky Way is believed to have four. The gas component is not included in the simulation.
To calculate the matter density distribution of the Galaxy (which we assume to be the DWD number density distribution) we use the kernel density estimation (KDE) method. We approximate the mass-density distribution of the Galaxy by placing a function called a kernel, in our case a three dimensional Gaussian, at the position of every simulation particle. The superposition of these functions is then used as a probability density function to populate the simulated Galaxy with DWD binaries. In other words, this is equivalent to assigning a DWD's position around a random simulation particle with a Gaussian probability density centred on the particle. We use a bandwidth of the KDE (equivalently the standard deviation of the Gaussian) of 10 pc. This choice minimises the residuals between the original simulation particle distribution and that reconstructed (backwards) from the KDE matter density distribution. In addition, this bandwidth allows us realistically describe the white dwarf's spatial distribution, whose local space density is 0.0045 pc −3 (e.g., Hollands et al. 2018).

GW emission from simulated binaries
Binary population synthesis provide us orbital periods P orb and component's masses m 1 and m 2 , while the Galaxy simulation gives us the DWDs' 3D positions. Here we provide relations between these quantities and parameters characterising GW signals of DWD binaries: GW frequency f , the time derivative of GW frequency f , GW amplitude A, sky position in LISA's reference frame (θ,φ), inclination ι, polarisation angle ψ, and initial orbital phase φ 0 . The first three are obtained from binary properties as is the triad of Cartesian coordinates centred on the Galactic centre. Included are the bulge and disc components.
where M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 is the chirp mass and d is the luminosity distance, G and c are respectively the gravitational constant and the speed of light. We compute sky coordinates and the distance to the source in the LISA reference frame as with (x ecl , y ecl , z ecl ) being the triad of Cartesian coordinates in the heliocentric ecliptic system and d = x 2 ecl + y 2 ecl + z 2 ecl . We assume the distance of the Sun from the Galactic centre is 8.1 kpc (Gravity Collaboration et al. 2019;Reid et al. 2019) and a viewing angle with respect to the bar's long-axis is 30 • (Cao et al. 2013). Finally, binary inclination angles are drawn from a uniform distribution in cos ι, while ψ and φ 0 are randomly drawn from uniform distributions respectively between [0, π] and [0, 2π].

MOCK GW OBSERVATIONS
In this work we consider the latest LISA mission concept consisting of three identical spacecrafts in an equilateral triangle configuration of 2.5×10 6 km per side (Amaro-Seoane et al. 2017). The spacecrafts are designed to exchange laser links (2 links per arm of the triangle) closing the triangular configurations. This design allows generation of two sets of data streams yielding two independent time-series with uncorrelated noise, in this way maximising the signal-to-noise ratio (SNR) of a GW event (Vallisneri 2005). We assume the nominal detection threshold corresponding to a signal-to-noise ratio (SNR) of 7, a mission duration of 4 years and instrument noise model from the LISA mission proposal Amaro-Seoane et al. (2017).
To compute the number of LISA detections and predict the uncertainties on their parameters we employ the Mock LISA Data Challenge (MLDC) pipeline developed by Littenberg et al. (2013) 1 and more recently updated by A. Petiteau. The pipeline is designed for optimising the analysis of a large number of quasimonochromatic GW sources simultaneously present in the data. To extract LISA detections from the input catalogues the pipeline iteratively computes the overall noise level by running a median smoothing function through the power spectrum of the population and extracts high SNR sources until the convergence. Finally, the pipeline computes the errors on GW observables by computing Fisher information matrices for detected sources.
We find 21.7 × 10 3 DWDs with SNR > 7, i.e. individually resolved by LISA, in agreement with our previous results (Korol et al. 2019). We show their number density maps projected in the Galactic plane in Fig. 3; we indicate the position of the LISA detector with the black triangle.
In contrast to electromagnetic studies, GW observations provide measurements of the amplitude of the waves. Even though the amplitude scales as 1/d, shallower than the intensity 1/d 2 scaling, one still detects fewer sources at larger distances. This bias can be corrected assuming that the population of DWDs is homogeneous in its binary properties throughout the Galaxy and that the population evolves in frequency only due to GW emission. Under these assumptions the fraction of observed sources with respect to the total (input catalogue) becomes a function of the distance only. We follow the approach outlined in Korol et al. (2019, eqs. (B1)-(B3)) to construct a de-biasing function. First, we compute the detection fraction as a function of distance using the full DWD catalogue (detected and undetected), and we fit the resulting distribution with a power law. Then, the multiplicative inverse of this power law is used to assign a weight to each detected source. Effectively, each weight represents the number of non-detections for each detected source. The power law obtained from our data is: The effect of the correction is clearly visible in the bottom panel of Fig. 3, where there are relatively fewer sources around the Sun with respect to the upper panel.
The two maps in Fig. 3 reveal the potential of using GW sources as tracers of Galactic structure that allow the reconstruction of the From geometrical considerations we can see that φ c (R) = arcsin H 2R whole surface of the disc including the far side of the Milky Way. In particular, the denser central Galactic regions are well mapped and the shape of the bar is already clearly visible without further data processing, unlike for optical observations (e.g., Anders et al. 2019). On the other hand, just hints of the spiral structure are present at distances < 5 kpc from the Galactic Centre but overall it cannot be clearly identified from visual inspection alone.

GALAXY STRUCTURE ANALYSIS
In this section we outline the method adopted to recover the properties of our fiducial Galaxy model from mock LISA observations. We carry out the analysis of the density distribution in Fourier space. First, we give the definition of the Fourier modes and discuss their meaning in our study. Next, we consider simple analytic distributions useful for understanding results presented in the next Section.
To define Fourier modes of the Milky Way's density distribution, we integrate over concentric rings centred on the Galactic Centre. Within these rings, both the bar and spiral arms are expected to appear as periodic over-densities. We define the Fourier modes A m as where m denotes the angular frequency of the mode, R is the cylindrical Galactocentric radius, φ is the polar angle in the Galactic plane, and n (R, φ) is the number density distribution projected on the Galactic plane. A prominent A 2 mode implies a pair of overdensities and thus implies either a bar or a pair of spiral arms. The difference between the bar and spiral arm structures is encoded in the mode (complex) phase. A bar has a complex phase which is constant with R, while spiral arms are characterised by a radial dependent phase. In the following we normalise the Fourier modes A 2 to the zero mode A 0 , representing the total mass in the ring, so that the ratio is confined between [0, 1]. Next, we analytically derive the Fourier modes for a number of simple analytic density distributions: a uniform density bar, logarithmic spirals, a toy barred spiral galaxy and a 2D Gaussian density bar. These models are useful for interpreting the results when analysing the density distribution of the DWDs detected by LISA.

Uniform density bar
Let us consider a bar with a uniform density n 0 , half-length R b and full width H. The density profile can be expressed as 6 Martijn J. C. Wilhelm Figure 5. Three examples of density distributions in left panels and respective the trend of the normalised m = 2 mode magnitude (solid lines) and (complex) phase (dotted lines) with R in right panels. Top panels: A uniform density bar under two different rotations. The magnitude vanishes for radii smaller than half the bar width, and then quickly converges to 1. The phase is constant for all radii, and is shifted by twice the angle to the x-axis. Middle panels: A toy barred spiral with infinitesimally thin bar and logarithmic spirals. In this case the magnitude is constant for all radii and equal to 1, while the phase is constant for radii smaller than the bar length and drops beyond. The jumps in the phase are due to the phase-wrapping to the interval (−π, π]. Bottom panels: A 2D Gaussian bar density model. The grey solid line denotes the profile of a bar with half the width of the model shown. The magnitude is 0 only at the origin and converges to 1. The speed of this convergence depends on the bar width. The phase is constant at 0 (except at the origin where it is undefined) like the uniform bar.
For certain radii R there is a critical angle φ c (R) = arcsin(H/2R) that indicates the boundaries of the bar in polar coordinates. This is represented in Fig. 4. Using the definition in Eq. (8) we can write the Fourier modes as This only holds for m > 0, and the absolute value of the radial coordinate H/2 < |R| < R b . Inside H/2 and far outside R b the Fourier modes vanish. Note that there is an intermediate regime between R = R b , where the circle with radius R covers the bar from −φ c to φ c , and R = R 2 b + (H/2) 2 , where the bar is entirely inside the circle (and the Fourier modes vanish). We note that in this regime the integration boundaries will have a more complex expression; although a detailed description is beyond the scope of the paper.
The m = 0 mode is so the normalised Fourier m = 2 mode is Rotation of the bar by an angle φ 0 implies a shift in the integration boundaries by the same angle. This results in a constant (complex) phase factor for the Fourier modes, i.e.
See the top left panel of Fig. 5 to visualise the density profile of a uniform density bar with φ 0 = 0 in blue and φ = π/4 in red. The top right panel of Fig. 5 shows the normalised m = 2 Fourier mode (solid black line) and the phase (dotted lines) for two configurations. For R < H/2 the magnitude is zero and for R > H/2 it quickly rises and converges to 1 at R = R b . In contrast, the phase of the mode remains constant at all radii with a value equal to 2φ 0 as discussed above.

Logarithmic spirals
Although both simulations and observations do not exactly follow logarithmic spiral structure, here we consider this example as a helpful analytic tool to construct a simple barred spiral model in Section 4.3. The logarithmic spiral can be parametrised as where R 0 is the distance to the origin for φ = 0 and fixes the phase of the spiral, and b is related to the pitch angle. The pitch angle is defined as µ = arctan | 1 R dR dφ |, thus it can be expressed in terms of b as µ = arctan |b|. The density distribution of a set of N evenly distributed, infinitesimally thin logarithmic spirals is then Again, using Eq. (8) we can express the m = 2 mode of this density distribution as From this equation it follows that the normalised m = 2 Fourier mode is

Barred spiral
Now we combine the uniform density bar and the logarithmic spirals to construct a toy model for a barred spiral galaxy. We consider the bar to be infinitesimally thin (i.e. φ c → 0) out to some radius R b and a pair of logarithmic spirals starting from the edges of the bar (i.e. we set R 0 = R b ). The magnitude of the normalised m = 2 Fourier mode is 1 and its phase Φ is of the form From which follows that rotating this toy galaxy by an angle of φ 0 results in a phase shift of 2φ 0 (considering the bar example and continuity). This is represented in the middle panels of Fig. 5. The middle right panel shows that the magnitude of the normalised m = 2 mode (solid line) is equal to 1 at all radii, while the phase (dotted lines) is constant for R < R b and drops afterwards. The discontinuities in the phase are due to its wrapping to the interval (−π, π].

Gaussian bar
Finally, we consider a more realistic example of a bar with density described according to a 2D Gaussian profile. The Fourier modes for a Gaussian bar can we written as where k = (R x /R y ) 2 with R x and R y being characteristic scales in the Xand Y -directions. Useful for later, we define the bar length R b in this model as while the bar full width H is Consequently, the axis ratio H/2R b is 1/ √ k for k > 1, and √ k if k < 1. Note that a bar rotated by 90 • has R x and R y flipped, but it will have the same Fourier magnitude profile as the rotation is encoded in the phase.
The integral in Eq. (19) does not have an analytic solution and needs to be integrated numerically. We plot the numerical results in the lower panels of Fig. 5. The grey and the black solid lines in the right panel represent the trend of | A 2 /A 0 | with R respectively for the case of half bar width and full bar width. Both start from Phase (A 2 (R)) (rad) Figure 6. The magnitude and phase of the normalised Fourier mode | A 2 /A 0 | of both the total stellar population (blue line, labelled "total density profile") and the population of individually detected DWDs (black points, labelled "GW data"), as a function of Galactocentric radius. The horizontal, red dotted line denotes a phase of 60 • ; this is the phase expected for the m = 2 mode of a bar rotated by 30 • (see Section 4). The vertical, red dash-dotted line denotes a radius of 4.5 kpc, the fiducial bar length. 0 in the origin and converge to 1 at 2 characteristic scale radius (R x and R y ). It is evident from the two examples that the speed the convergence depends on the width of the bar. The phase (dotted line) is constant at 0 (except at the origin) as in the case of the uniform bar (top panels in Fig. 5).

RESULTS
We calculate the number density distribution in Galactocentric coordinates of gravitational wave sources (cfr. Fig. 3 lower panel) as a superposition of delta functions weighted by the de-biasing factor w at the location of each detection, This is the formal definition of the density distribution of a collection of point particles. The sum is performed over all detections i, and each is weighted by the de-biasing weight, which depends on its position through the distance from the detector. Entering this into In practise, we calculate Fourier modes (Eq. 23) by radially binning all detections to make a histogram, where every detection is weighted by its de-biasing weight and a complex phase factor. We choose a bin size [R, R + dR) with dR = 100 pc in order to obtain the mock GW data shown in Fig. 6 (black circles).
To quantify error bars, we then repeat the above procedure for 10 4 realisations of the binary positions in the Galaxy. Specifically, each realisation is obtained by randomly drawing distances and sky positions from Gaussian distributions centred on the measured values, with standard deviation equal to the measurement errors provided by the MLDC pipeline. This procedure allows us to assign statistical uncertainties to a given realisation, calculated as the standard deviation over all data realisations (error bars in Fig. 6). In particular, we show the normalised m = 2 Fourier mode's magnitude (top panel) and phase (bottom panel) as a function of the Galactocentric cylindrical radius. For comparison we also represent the Fourier transform of the total stellar distribution in the Galaxy with the blue solid line. Note that the blue line in the top panel shows two peaks: one around 2.5 kpc and another one around 6.5 kpc. The first peak is caused by the central bar, the second one by the spiral arms. The mock data follow the total number density profile well up to roughly 7 kpc, describing well the first peak. Beyond 7 kpc the data become too noisy to detect over-densities.
The bottom panel of Fig. 6 shows the phase of the total density profile (blue line). It stays at a constant value of π/3 up to ∼ 4 kpc, before decreasing down to −π at ∼ 7 kpc, where it suddenly jumps to π and then it decreases again. As shown in the middle panels of Fig. 5, this behaviour is characteristic of the barred spiral galaxy. In particular, our model example in Sect. 4.3 reveals that the phase of the A 2 mode remains constant for the extent of the bar, while the jumps are due to the phase wrapping. The mock data in Fig. 6 closely follow the blue line up to roughly 5 kpc and become more scattered at greater radii. As can be seen in Fig. 3, this is due to the low density contrast between DWDs in the spiral arms and the background disc. Therefore, from visual inspection of Fig. 6 we can conclude that by using the reconstructed DWD density maps from LISA observations one can clearly identify the bar only. The spiral arms' identification and characterisation is more challenging than for the bar, whether we exploit the m = 2 Fourier magnitude or the phase.
To derive the bar's parameter, we fit the mock data with the Fourier transform of the 2D-Gaussian profile defined in Eq. (19). Specifically, we fit for the scale length R x and the square inverse axis ratio k. This choice is motivated by the fact that when computing the Fourier magnitude of the Gaussian profile, the two scale lengths are degenerate: switching the scale lengths -corresponding to a rotation of 90 • -does not influence the Fourier magnitude, as mentioned in Sect. 4. The fit is performed from the origin up to and including the maximum magnitude of A 2 /A 0 . In both cases k > 1, meaning that the bar length is computed from the bar width and axis ratio, which means that its error is the compound error of those of the other two values, in turn leading to a larger error.
We then recover an axis ratio of 0.27 ± 0.03, R b = 4.34 ± 0.51 kpc, and H = 2.32 ± 0.10 kpc. A similar fit to the Fourier magnitude profile of the stellar mass distribution yields an axis ratio of 0.28 ± 0.03, R b = 4.53 ± 0.59 kpc, and H = 2.52 ± 0.16 kpc. Overall, errors vary from 5 -15 per cent and these 'observed' and 'true' values for the bar axis ratio and length agree well within 1σ. Instead, the observed value for the bar width H differs by about 2σ with respect to the true value. Although we underestimate both the bar length and bar width, the axis ratio remains consistent with the true value well within 1σ.
As discussed in Sec. 4, bar length and orientation can be independently recovered from the phase of the m = 2 mode. This is possible using the fact that the phase stays constant at Φ (A 2 ) = 2φ 0 for the extent of the bar and deviates from this value beyond it (see Fig. 6). To perform these measurements we apply the following criterion. Between 0.5 kpc (to exclude the bulge) and R we fit a horizontal line Φ (A 2 ) = const and compute the χ 2 value of the fit up to that given point. In this way, we obtain a function χ 2 (R) that quantifies the goodness of fit to a constant value, for data up to R. Using the statistical distribution of χ 2 we can then compute the cumulative distribution function (CDF) for each fit, that can be interpreted as the probability that a randomly generated set of data fits the model better. We compute 1 -CDF, i.e. the probability of generating a random set of data that fits the model worse. We then have a function that quantifies, for every radius R, how well a constant can fit points up to R. The bar length is then defined as the R up to which the fit is good. We take a lower limit of 0.05 2 for an acceptable fit. We find that for R 4.25 kpc, our statistic is always greater than 0.6, while for R > 4.25 kpc, it is always smaller than 10 −6 . The boundary thus appears well-defined, and relatively insensitive to our choice of the threshold probability value. In this way we recover a bar length of R b = 4.25 kpc. We retrieve the viewing angle by taking the best fit to the data up to the bar length, which gives us a value of φ 0 = 30.70 ± 0.83 • , in agreement with the true viewing angle of 30 • .
In Fig. 8, we visualise our results by over-plotting our bar width 2 i.e. the probability that a set of data randomly generated from the model fits the model worse than this data is 5 per cent (white dotted line) and the two bar length estimates (white and black dotted lines) onto the total density distribution in the Galactic plane. The shaded regions show the 1σ confidence level of our fit. The values for the bar length and width fully enclose the bar. We note that the two values of the bar length have been independently derived and the fact that they agree well with each other and with the true value within one sigma confirm the robustness of our method. To elaborate on the independence of the two measurements, the bar length recovered from the Fourier magnitude is a scale length, describing the decay of a density profile, while that recovered from the Fourier phase can be interpreted as a physical length, namely the extent of the feature that is constant in Fourier phase (and which feature we identify as the bar). The figure also shows the agreement between the actual viewing angle (yellow line) and our viewing angle measurement (black dashed line).
In order to investigate the sensitivity of our results to different viewing angles we perform a set of 13 simulations for a number of vantage points from 0 to 180 • . In all simulations we place the virtual LISA detector at 8.1 kpc from the Galactic Centre and rotate it around a half-circle in steps of 15 • . The top panel of figure 9 shows the results for the bar lengths, bar widths, and axis ratio's derived from the magnitude of the m = 2 mode for the different vantage points. The bar lengths and axis ratios are generally consistent within 1σ with the true value (red line) obtained from fitting to the full stellar density profile. The derived bar widths, instead, underestimate the true value by more than 1σ for most viewing angles. Likewise, we investigate the sensitivity to the viewing angle of the Fourier phase modelling (Fig. 9, bottom panel). The bar lengths show a discretisation due to our definition, which was based on choosing a discrete bin. All values obtained fall in one of three bins, meaning that they all fall within a range of 0.2 kpc. The fitted viewing angles are for the most part consistent with their true values, with errors on the order of 1 • .
For each viewing angle, the de-biasing function was recomputed, and we compare the functional parameters for each in Table 1. We note that their standard deviations are on the order of a few per cent and that their means and medians agree to about 10 per cent of a standard deviation, implying that these are robust results. These values can potentially be used for as de-biasing function for actual LISA observations of DWDs. However, the degree to which our debiasing function approximates the true one depends on the quality of the LISA noise model and on the DWD population model, and warrants a separate more careful investigation.

DISCUSSION AND CONCLUSIONS
It is well established that the Milky Way presents a complex baryonic structure including bulge, a thin and a thick disc and a stellar halo (see Bland-Hawthorn & Gerhard 2016, for a review). In particular, the structure of the inner Galaxy is not well known due to heavy extinction and stellar crowding. The best available structural information on the stellar population of the inner Galaxy comes from large samples of Red Clump Giant stars visible in near-infrared band, for which individual distances can be determined (Minniti et al. 2010;Wegg & Gerhard 2013;Simion et al. 2017). Numerical simulations show that bulges represent the inner part of a longer, planar bar structure that formed through buckling out of the galaxy plane and/or from orbits in vertical resonance (e.g., Athanassoula 2005). Therefore our Galaxy is also expected to have a thin bar component extending well outside the bulge. However, finding the Galactic planar bar and characterising its properties has proven difficult, again because of intervening dust extinction and the superposition with the star-forming disc at low latitudes toward the inner Galaxy.
In this paper, we exploit the fact that GW observations do not suffer from these effects and that GW strength decreases slower with distance than light, allowing us to reach the far side of the Galaxy. In particular, we propose and develop for the first time a methodology for investigating the inner structure of the Milky Way using future LISA observations. We combine binary population synthesis techniques with a high resolution Milky Way simulation for a more realistic modelling of the Galactic structure (D'Onghia & Aguerri 2019).
Our results reveal that the shape of the bar is already apparent by simply plotting the 3D positions of LISA detections (unlike with optical observations Anders et al. (2019)), while the spiral arms remain elusive. To recover the bar's parameters, we analyse the density distribution in Fourier space of LISA detections. We find that the m = 2 Fourier mode is the strongest, implying the presence of a bar plus a pair of spiral arms structure.
By fitting a 2D Gaussian density function to the Fourier transform of the bar, we recover a length R b = 4.34 ± 0.51 kpc, a width H = 2.32±0.10 kpc and an axis ratio 0.27±0.03. The results for the bar length and bar axis ratio are consistent within 1σ with the true values, obtained by performing the same analysis on the total stellar population. The bar width, instead, is underestimated by ∼ 2σ, although the best-fit and true values differ by just ∼ 8 per cent. We note that the bar length was also underestimated, though within 1σ, and that the axis ratio agrees to within 1σ and 4 per cent of the true value. Bar parameters obtained assuming different viewing angles are consistent with each other within 1σ, and therefore similarly compare to the true values. By fitting the phase of the m = 2 mode as a function of R, we recover a viewing angle φ 0 = 30.70 ± 0.83 • , consistent with the true angle within 1σ. Viewing angle estimates inferred from other vantage points are generally consistent with the true value within 1σ, with the fitting error being on the order of a degree. Likewise, the bar lengths recovered with this method for different viewing angles do not vary by more than 4 per cent.
Let's now compare our results with the state of the art for these measurements. Electromagnetic studies of the geometric properties of the Galactic bar have used the optical and infra-red stellar tracers. Recent examples include infrared observations of stars using the 2MASS survey (Robin et al. 2012), optical observations of Red clump stars (Cao et al. 2013) and RR Lyrae stars (Pietrukowicz et al. 2015) from the OGLE survey, and optical observations of red clump stars with the VVV survey (Wegg & Gerhard 2013;Simion et al. 2017). The geometric parameters of the bar obtained by these five studies vary between 10 • − 45 • and 0.25 − 0.6 in viewing angle and in-plane axis ratio respectively (Simion et al. 2017, figure 17). Error estimates are reported to be of the order of a few degrees for the viewing angle (Robin et al. 2012;Pietrukowicz et al. 2015;Wegg & Gerhard 2013); around tens of per cent for scale lengths (Robin et al. 2012); and around~5 per cent for the in-plane axis ratio (Pietrukowicz et al. 2015). Instead from GWs measurements, one can recover the bar inplane axis-ratio, scale length and viewing angle accurately (within 1σ from the true values) and with a precision of < 10 per cent for the two former quantities, and 1 • for the latter one. Therefore, GWs promise to be a competitive observational window for the central region of the Galaxy.
Nevertheless, the analysis in this paper should be refined in a few ways. We have assumed a simple Galactic plane-projected, bi-axial Gaussian profile, whereas the bar of the Milky Way is more peanut-shaped and with a vertical structure we have currently ignored. Additionally, we didn't attempt to simultaneously recover the scale height of the bulge and disc, although possible as shown by Benacquista & Holley-Bockelmann (2006); Korol et al. (2019); Breivik et al. (2019b), albeit with simple analytical density models. Another improvement pertains to our choice for the star formation rate. While we assumed that the bar follows the same star formation as the disc, one should instead self-consistently simulate the star formation history of the bar itself. The different star formation histories between the disc and the bulge and bar can potentially give rise to two distinct populations. Finally, the future is currently offering at least two more space-based interferometers beside LISA: the DECI-hertz inteferometer Gravitational wave Observatory (DE-CIGO, Sato et al. 2017) and TianQin (Luo et al. 2016). An analysis that would consider the combined detections would likely result in a more effective and precise tomography of the Galaxy.
Despite these limitations, our analysis represents a first proof of concept of the feasibility to recover and study the inner parts of our Galaxy using GWs, and we anticipate that a synergistic study with more traditional electromagnetic tracers will be a new and very effective Galactic investigation strategy in the 2030s.