Hic sunt dracones: Cartography of the Milky Way spiral arms and bar resonances with Gaia Data Release 2

In this paper we introduce a new method for analysing Milky Way phase-space which allows us to reveal the imprint left by the Milky Way bar and spiral arms on the stars with full phase-space data in Gaia Data Release 2. The unprecedented quality and extended spatial coverage of these data allowed us to discover six prominent stellar density structures in the disc to a distance of 5 kpc from the Sun. Four of these structures correspond to the spiral arms detected previously in the gas and young stars (Scutum-Centaurus, Sagittarius, Local, and Perseus). The remaining two are associated with the main resonances of the Milky Way bar where corotation is placed at around 6.2 kpc and the outer Lindblad resonance beyond the solar radius, at around 9 kpc. For the first time we provide evidence of the imprint left by spiral arms and resonances in the stellar densities not relying on a specific tracer, through enhancing the signatures left by these asymmetries. Our method offers new avenues for studying how the stellar populations in our Galaxy are shaped.


Introduction
A long-standing problem in galactic dynamics is how gravitational instabilities of stellar discs shape the large-scale morphology of galaxies, and in particular their conspicuous spiral arms. While the spectacular images of spiral galaxies mainly reflect young stars formed from interstellar gas in spiral arms, all disc stars can cross the arms and contribute to their overdensities; therefore, the entire disc should bear the imprint of its spiral structure. In the Milky Way, the position of our Solar System near the plane of our Galaxy, distance uncertainties, and interstellar extinction have complicated attempts to learn directly about the large-scale disc morphology. Using observations of massive stars and pulsars (Taylor & Cordes 1993;Xu et al. 2006Xu et al. , 2018, HII regions (Georgelin & Georgelin 1976;Russeil 2003), and masers (Reid et al. 2014(Reid et al. , 2019Sanna et al. 2017) numerous studies have attempted to delineate the spiral arms of our Galaxy (Vallée 2014). However, the structure of the disc in old stars, beyond the solar neighbourhood, remains largely terra incognita.
The second data release (DR2) of the ESA Gaia mission has provided the largest available 6D phase-space (positions and velocities) dataset for 7.2 million stars brighter than G RVS = 12 mag (Gaia Collaboration et al. 2018b), making possible precise studies of the Milky Way structure and kinematics on large scales. Gaia data have already revealed signatures of nonequilibrium and ongoing vertical phase mixing in the Milky Way disc , likely induced by a previous pericentric passage of the Sagittarius dwarf galaxy  Here be dragons, a phrase famous in medieval cartography when dragons and sea monsters were used to designate uncharted and possibly dangerous regions. Bland-Hawthorn et al. 2019;Haines et al. 2019) or by internally driven bending waves (Khoperskov et al. 2019;Darling & Widrow 2019). A large number of kinematic arches with various morphologies not known prior to Gaia were found (Ramos et al. 2018;Kushniruk & Bensby 2019;Monari et al. 2019) and largescale wiggles were discovered in the V φ − R plane Antoja et al. 2018). These features can be interpreted as the signature of the impact of an external perturbation (D'Onghia et al. 2016;Laporte et al. 2019) or as evidence of the spiral structure and the bar of the Milky Way (Hunt et al. 2018;Fragkoudi et al. 2019). At the same time, in the solar vicinity the Galactic spiral structure is believed to generate various patterns in the phase-space (Siebert et al. 2012;Williams et al. 2013), but uncertainties on the location and strength of the spiral arms have so far prevented us from explaining the observed patterns. In this work, we use the high-quality Gaia DR2 sample of stars with radial velocities (hereafter GRV2), together with a new method to highlight the stellar density structures, to explore the Milky Way spiral arms without relying on any specific stellar tracer. The outline of the Letter is as follows: in Sect. 2 we describe the data and the method we used; in Sect. 3 we examine the phase-space structure of the Milky Way and compare its various features with Milky Way-type galaxy simulations; and in Sect. 4 we summarize the conclusions that can be drawn from our study. A&A 634, L8 (2020) 2018a). From this sample we selected stars with positive parallaxes, and relative errors on parallaxes of less than 10%. Distances were computed by inverting parallaxes. For calculating positions and velocities in the galactocentric rest-frame, we assumed an in-plane distance of the Sun from the Galactic centre of 8.19 kpc (Gravity Collaboration 2018), a velocity of the Local Standard of Rest, V LSR = 240 km s −1 (Reid et al. 2014), and a peculiar velocity of the Sun with respect to the LSR, U = 11.1 km s −1 , V = 12.24 km s −1 , W = 7.25 km s −1 (Schönrich et al. 2010). After applying these selections, our final sample consists of ≈5.3 M stars, whose density distribution, projected onto the Galactic plane, is shown in Fig. 1a. As expected, the main shape of the distribution is controlled by the Gaia selection function, with the density rapidly decreasing with the distance from the Sun. We also investigated the sensitivity of our results on the Gaia parallax biases (Schönrich et al. 2019) and did not find any systematic changes in the structures we discovered.
The detection of stellar structures in the Galaxy faces two problems of different natures. On the one hand, the epicyclic motion of stars, which is the radial oscillation around their guiding centres, blurs all kinds of morphological features in discs. This effect is clearly seen in galactic structures made of multiple stellar populations with different kinematics, making any density structure more evident for younger, dynamically colder stellar populations with lower orbital eccentricities (Debattista et al. 2017;Khoperskov et al. 2018). On the other hand, the Gaia and Gaia RVS limiting magnitudes, the shape of the Gaia footprint, and the effect of extinction have the consequence of limiting the accessible portion of the disc to a few kpc from the Sun, with a density of sources that decreases rapidly with distance from the Sun (see Fig. 1a). To tackle the first of these two problems, it is convenient to analyse stellar distribution in what we call the "guiding coordinate space", i.e.
where R g = L z /V LSR is the guiding radius normalized to the LSR velocity taken to be V LSR = 240 km s −1 , L z is the instantaneous angular momentum of the star L z = R × V φ , R its galactocentric distance, φ is the azimuthal angle around the galactic centre clockwise from the direction towards the Sun, and V φ its azimuthal velocity in the Galactic plane. By transforming to (X g , Y g ) or (R g , φ) phase-space, we change the radial positions of stars without changing their azimuths φ. This allows us to minimize the effect of the radial oscillations of stars around their guiding centres, and hence sharpens the face-on features in the stellar density distribution. The guiding space coordinate transformation is based on the approximate conservation of angular momentum along the normal to the disc. This assumption is supported by the fact that the velocity field in the region from outside the bar to beyond the Sun is not far from circular (Reid et al. 2014;Bovy et al. 2015;Bland-Hawthorn & Gerhard 2016), despite perturbations from the bar and spiral arms, and the vertical disturbances of the disc discovered recently. In order to reduce the effect of the Gaia DR2 selection function, we apply an approach based on a re-sampling of the faceon density map in the spatial coordinates. First, we introduce a uniform Cartesian grid in spatial coordinates (X, Y) with the cell size of 100 × 100 pc 2 . Then we populate each grid cell by 100 randomly selected stars located within this cell. Such a single sampling procedure results in a constant number of stars in the (X, Y)-plane on the Cartesian grid. Then by repeating random re-sampling for each cell of the (X, Y)-grid we generate 100 independent selections of stars from the entire Gaia DR2, where each subsample contains 100 × 100 × 100 = 10 6 stars. Thanks to the large statistics of the Gaia DR2 catalogue, we can repopulate our (X, Y)-grid by different subsamples of stars, with repetitions, especially far beyond the solar vicinity where the number of the observed stars is lower. Then each subsample of stars (with a flat distribution in (X, Y)-coordinates) is converted into (X g , Y g ) guiding phase-space map (see Eq. (1)) where the flat distribution in (X, Y) plane transforms into a map depicting several large-scale features.

Results
Figure 1 (top) shows the mean density map ρ(X g , Y g ) for 100 re-sampled selections of stars where several horizontally extended overdensity features can be found across the entire disc area covered by Gaia DR2. To sharpen even more the features seen in the ρ(X g , Y g ) map, we produce a residual map by using an unsharp masking approach. In Fig. 1 (bottom), we show the residual map defined as where ρ(X g , Y g ) denotes the mean number density map constructed by convolving the density distribution (X g , Y g ) with a 2D Gaussian kernel of width 400 pc. The combined re-sampling technique and transformation to guiding coordinates reveal several new stellar density structures. Figure 1b shows the stellar density now plotted in the guiding coordinate space (X g , Y g ). Large-scale overdensities of length 3−8 kpc in guiding coordinates cover most of the disc observed by Gaia. Structures are seen even more clearly once the unsharp masking is applied. The result is presented in Fig. 1d where six outstanding narrow overdensities of ≈5−12% amplitude are visible. The most prominent features are seen close to the solar radius: one crosses the line from the Sun to the Galactic centre at Y g ≈ −7 kpc, and another passes through the solar vicinity at Y g ≈ −8 kpc. Two additional structures are detected in the inner galaxy, between Y g = −4.5 and Y g = −6 kpc. In the outer disc, we find an extended structure between Y g = −9 and Y g = −10.5 kpc and a smaller one just beyond the solar radius between Y g = −8 to Y g = −9 kpc. Coloured lines in Figs. 1c,e highlight the main structures that appear in the guiding coordinate space.
To test the possible effects of extinction, we investigated a cross-matched sample between GRV2 and WISE (Wilson & Naylor 2018). We selected stars in this sample with absolute W1 magnitude in [−1,−2], and considered the redder and bluer half in G-W1 colour separately. We found that the redder half of this sample, expected to contain a larger fraction of reddened stars, showed a significantly noisier representation of structures in guiding space than the bluer half, as well as a dominant extinction feature in the direction of the Galactic centre. This suggests that the guiding space structures seen in GRV2 are not the result of sample biases due to extinction, but that, on the contrary, extinction likely reduces the amplitudes of the structures we found. Therefore, our estimates of their amplitudes should be considered as lower limits.
In order to understand what causes these structures we consider two Milky Way-type galaxy simulations. The first is an Nbody/SPH simulation of a spiral galaxy (Model I, see details of the model in Appendix A) in which the spiral structure in the stellar disc is weak and blurred by the epicyclic motions of the stars (see Fig. 2a). However, transforming the positions of the simulated stars to guiding coordinate space using Eq. (1) allows us to recover narrow spiral arm-like features (see Fig. 2b) whose morphology is in much closer agreement with the gas surface density map than the original stellar density (see Appendix A). For the L8, page 2 of 9 SDS1 SDS2 SDS4 SDS5 SDS6 SDS3 Fig. 1. Face-on view on the Milky Way disc with Gaia DR2: spatial and guiding coordinate space. Panel a: two-dimensional map in the Galactic plane of the number of Gaia DR2 stars with radial velocities, colour-coded by number of stars. The grey oval represents the orientation of the Milky Way bar (Wegg et al. 2015), and the solar position is shown by a yellow star. The stellar density distribution is shaped by the Gaia footprint and is heavily influenced by the fading of stars at greater distances and by the high interstellar extinction in the Galactic plane. Panels b and c: twodimensional map in guiding coordinate space of the mean number of stars over 100 random samples drawn from the homogenized spatial density map ρ(X g , Y g ). Panels d and e: unsharp masking map (see Eq. (2)) of the guiding space showing six newly discovered overdensities prominent in a wide area of the Milky Way disc, each containing between 2000 and 6000 stars, about 5-12% of the respective background. This estimate is a lower limit for the amplitude of the overdensities because of the dust extinction in the disc plane blurring their visible structure. These features are indicated in panels c and e, and are interpreted below as the imprint of the Galactic bar (corotation and outer Lindblad resonances) and the four-armed spiral structure in the Milky Way. The irregular vertical features near X ≈ 0 line are artefacts of the Gaia DR2 footprint in this area. simulated stars in each of these features we calculate the distribution of radial velocities (see Fig. 2c; more details in Appendix A). In each histogram we see a systematic inward motion, with the radial velocity distribution skewed towards negative values. This pattern is in agreement with predictions of density wave theory where, within corotation, the stars in the spiral arms move inwards towards the galactic centre (Lin et al. 1969;Siebert et al. 2012;Faure et al. 2014;Pasetto et al. 2016;Monari et al. 2016).
Another non-axisymmetric structure that could generate features in the guiding coordinate space is a stellar bar which gives rise to a number of resonances in the outer disc. In Figs. 2d-f we show the results from a snapshot of an N-body simulation of a barred galaxy without prominent spiral structure (Model II; see details in Appendix A). Here, two narrow quasi-circular overdensities appear in the guiding coordinate space outside the bar which are absent in the spatial density map. Spectral analysis of the orbits of particles in these overdensities confirms that these features are at the bar's corotation and outer Lindblad (OLR) resonances (see Appendix A for details). The observed kinematics of stars at the bar resonances depend on both the bar orientation and the distribution of stellar orbits. In particular, some stars near corotation follow orbits elongated along the bar major axis, showing periodic variations of their radial velocities. The bar OLR is populated by two orbital families with proportions varying along the OLR region. When the bar is oriented as it is in the Milky Way, and an area of the simulated disc is considered similar to that available with GRV2 (Y < 0, |X| < 5 kpc), we find that the OLR region is dominated by x 1 (2) orbits (Fragkoudi et al. 2019) that move outwards. This makes the mean radial velocity as seen from the Sun positive, with a velocity distribution skewed towards positive values. At the same time stars near corotation show different kinematics with more symmetric radial velocity distribution (see Fig. 2f and Appendix A for details).
We conclude from these theoretical models that both the spiral arms and major bar resonances are imprinted in velocity space. The extended spatial coverage, high sampling (even far from the solar vicinity), and unprecedented precision of the Gaia DR2 data allow us to see these signatures in the observed velocities. Figure 2g shows that four of the overdensity struc-tures obtained in Fig. 1 with the guiding space analysis (SDS 1, SDS 3, SDS 4, and SDS 6) have systematically negative radial velocity with a radial velocity distribution skewed towards negative values, just as the spiral arms in the simulated galaxy of Model I. Meanwhile SDS 2 has a wide and more symmetric distribution with only a small mean negative velocity, and stars in SDS 5 tend to move outwards as seen from the Sun, with the velocity distribution skewed towards positive values (see Fig. 2h), very similar to the barred galaxy simulation Model II.
We can also transfer the structures in Fig. 1 back to original Galactic plane coordinates (see Appendix B). To place all structures in the Galactic plane, we estimate their mean location by using Gaussian fitting along the galactocentric radial direction with 1D Gaussian models. Figure 3a shows the inferred locations of the Milky Way spiral arms (SDS 1,3,4,6). Overplotted are the positions of known high-mass star-forming regions (HMS-FRs) with accurate VLBI maser distances (Reid et al. 2014;Xu et al. 2016) which trace the Milky Way spiral arms in gas and young stars. These sample the Perseus arm (purple), the Local arm (orange), the Sagittarius arm (blue), and the Scutum-Centaurus arm (turquoise). The location of these HMSFRs is in excellent agreement with the stellar overdensities we find. The Local arm and the Perseus arm perfectly overlap with SDS 4 and SDS 6, the fourth and the sixth (from the Galactic centre) stellar density structures discovered in this study. The star-forming regions in the inner disc (associated with the Sagittarius arm and the Scutum-Centaurus arm) are spatially coincident with SDS 1 and SDS 3, the first and third of the discovered stellar structures; a few HMSFRs connecting SDS 1 and SDS 3 also overlap with SDS 2.
Fitting a logarithmic spiral arm model to the data in Fig. 3a, we measure the pitch angles (Davis et al. 2012) for the spiral arm structures (SDS 1, 3, 4, 6), from the inner to outer arm 14.7±1.1 • , 10.8 ± 0.8 • , 5.8 ± 0.7 • , and 7.7 ± 0.7 • , respectively. Therefore, our results suggest that the Milky Way has a four-armed spiral, in agreement with a number of studies of Milky Way spiral structure based on different tracers and techniques ( where the mean radial velocity within the arms is negative (inwards), exactly as expected from density wave theory (Lin et al. 1969). Panels d-f: Milky Way-type barred galaxy (Model II). Transforming stellar surface density (e) to guiding coordinate space (f) shows prominent overdensities at the bar's corotation (CR) and outer Lindblad (OLR) resonances, as determined from spectral analysis in Appendix A. Panel g shows a broad, nearly symmetric radial velocity distribution for stars at CR (blue), as viewed from the solar position at (x, y) = (0, −8.2) kpc, and an outwardskewed distribution at OLR (red). Panels g and h: kinematics of stars in Gaia GRV2 associated with the large-scale stellar overdensities (SDS) detected in guiding coordinate space (see Fig. 1). Panel g: radial velocity distribution for four structures (SDS 1, 3, 4, 6) reveals a mean radial motion towards the galactic centre similar to that seen in the spiral arms of Model I. Panel h: radial velocity distributions in SDS 2, 5 resemble the distributions at CR and OLR in the in the barred galaxy Model II. Error bars on all profiles are determined as described in Appendix A. Interestingly, the Local arm, which is the closest spiral arm to the Sun, appears to be a tightly wound structure, located between the two primary arms of the Milky Way, Sagittarius and Perseus. This result supports the idea that the Local arm is an intrinsic spiral arm, rather than an inter-arm spur as has been suggested for a long time (Elmegreen 1980;Scoville et al. 2001;La Vigne et al. 2006;Miyachi et al. 2019).
The SDS structures associated with the Milky Way bar resonances (SDS 2 and SDS 5) are shown as grey-shaded regions in Fig. 3a where corotation is placed at R ≈ 6.2 kpc and the OLR is just beyond the solar radius, at R ≈ 9.3 kpc. We note that the derived location of the resonances could correspond to a fraction of the entire resonance structure which is supported by our Model II, where both CR and OLR regions cover wide radially extended areas across the disc (see Figs. A.2 and B.1). This issue of completeness of the detected OLR stars in the disc beyond the Sun will be resolved by applying the approach presented here to stars with radial velocities in the forthcoming Gaia DR3. Although the Milky Way bar itself is not covered by our GRV2 sample, the location of the bar resonances allows us to estimate its angular rotation speed. For the significantly declining rotation curve of the Milky Way (Eilers et al. 2019) we obtain L8, page 4 of 9 values of 37−42 km s −1 kpc −1 in good agreement with recent in situ estimates from Milky Way bar kinematics and dynamics (Portail et al. 2017;Sanders et al. 2019;Bovy et al. 2019).

Conclusions
Understanding the global structure of our Milky Way, in all its constituents, is a classical problem of Galactic astronomy. Mapping the main non-axisymmetric features in the disc, as traced by their stars, is an important advance made possible only by the superb quality of the Gaia DR2 data; to date the Galactic spiral arms had only been detectable with young stellar populations and gas. In this paper we provide a means for identifying stars of spiral arms and bar resonances without relying on a specific stellar tracer, opening up several further avenues towards understanding the Milky Way.
For the first time we revealed the imprint left by the Milky Way bar and spiral density waves on the entire disc population; the Gaia Data Release 2 enabled us to recover the prominent fourarm spiral structure of the Galaxy (Scutum-Centaurus, Sagittarius, Local and Perseus spiral arms) to a distance of 5 kpc from the Sun (see Figs. 3b,c). We also independently constrained the Milky Way bar parameters by detecting the prominent signatures of its main resonances, which together with the properties of spiral arms provide us a much more complete picture of the asymmetries present in the disc of our Galaxy, and never seen before by Gaia. The presence of the bar and extended spiral structure evidently induce in-plane phase-space perturbations that imply that the Milky Way stellar disc is out of equilibrium on large scales, and most of the features in the velocity space are related to the presence of these non-axisymmetric structures in the disc.
In this work we also described a new idea for analysing the Milky Way stellar populations, which in combination with Gaia data and spectroscopic information (WEAVE, MOONS, 4MOST) will allow us to study stellar dynamics in the vicinity of the Galactic spiral arms, including radial migration of stars, as well as the influence of the Galactic bar on the disc, and thus the disc's dynamical evolution. It will also shed light on the chemical composition of stars in spiral arms and the role of these structures in shaping the stellar populations over time.
Acknowledgements. The authors thank Michael Fall for useful discussions during preparation of the work. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. PDM and MH thank the ANR (Agence Nationale de la Recherche) for its financial support through the MOD4Gaia project (ANR-15-CE31-0007, P.I.: P. Di Matteo). Development of the code for simulation of stellar-gaseous galactic discs by AK and SK was supported by the Ministry of Science and Higher Education of the Russian Federation (government task No. 2.852.2017/4.6). Numerical simulations were computed with support of the Russian Science Foundation (project no. 19-72-20089) by using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University (project RFMEFI62117X001). This work was also granted access to the HPC resources of CINES under the allocation 2017-040507 (Pi : P. Di Matteo) made by GENCI.

A.1. Model I: spiral galaxy model
In order to explore the nature of the large-scale overdensities we discovered in the guiding space of the Milky Way, we present an analysis of two different models of the Milky Way-type galaxies. In particular, to disentangle the perturbations of the guiding phase-space caused by the spiral arms and by the bar, we analyse a model of a multi-arm spiral galaxy (Model I) and a model with a bar and lack of spiral arms (Model II). Model I is an N-body/hydrodynamical simulation of the stellar-gaseous galactic disc of the Milky Way-type galaxy. The gaseous subsystem is modelled as Smoothed Particle Hydrodynamics (SPH), and the dynamics of the stellar disc is modelled by using a direct particle-particle integration scheme. To avoid the change in the bar size due to possible angular momentum exchange between the disc and the halo-bulge system resulting in a slowdown of the bar, we use a rigid galactic potential for spherical components. The importance of a live DM halo has been noted in several works, where the angular momentum exchange between the halo and the disc plays an important role in increasing the bar length over time (Debattista & Sellwood 2000;Athanassoula 2003). However, the impact of live DM halos on the formation and evolution of the spiral arms is not evident. Rigid halo simulations accurately reproduce the radial migration of stars (e.g. Grand et al. 2012Grand et al. , 2013Vera-Ciro et al. 2014) and are in agreement with theories of spiral structure formation (e.g. Sellwood 2011;Khoperskov et al. 2012a,b;D'Onghia et al. 2013;Sellwood & Carlberg 2019). Thereby we are confident that, at least on the qualitative level, our Model I is a reasonable guide for understanding the Milky Way spiral arms kinematics. In our simulation, the number of particles is 2 21 and 2 22 for stars and gas respectively, and the softening length is about 40 pc for both components. The generation of the model initial equilibrium follows an iterative method of the Jeans equation solution (Khoperskov et al. 2012a). The simulation was performed with a direct GPU N-body/SPH code which provides a very high accuracy of the numerical integration (Khrapov & Khoperskov 2017). We focus our analysis on a single snapshot of the evolution at 1.3 Gyr when an extended multi-arm tightly wound spiral structure formed in the model.
In Fig. A.1 we show this simulated multi-arm spiral galaxy where the spirals are mostly seen in the gas density distribution while the stellar density is smooth and does not reveal any clear signatures of the spirals arms. Then, as in the Gaia DR2 data analysis, we calculated the stellar density distribution in the guiding (X g , Y g ) phase-space coordinates where the density map (see Fig. A.1a) demonstrates a few large-scale features similar to those we discovered in Gaia DR2. To proceed further with simulated data, we selected stars along four major features in guiding phase-space and show the face on distribution of the stars associated with these features (see Fig. A.1e). These stars, taken in narrow regions in the guiding space, cover a large area in the disc plane, with some overlap, and in order to constrain the locations of the structures we provide a 1D Gaussian fitting of the structure locations and their width in (X, Y) coordinates, similar to those we use for the Gaia data. In Fig. A.1f we compare the derived stellar density structures location with the gas density distribution. We note that the stellar spiral arms mostly co-exist with the gaseous arms. Because the gas component does not have the epicyclic motion, which leads to a large variation in the velocities for stars, there is a narrow structure in the gas seen along the stellar arm. However, at the edge of spiral arms, a spatial offset between gaseous and stellar spirals is clearly visible, which is expected and is due to hydrodynamical instabilities of the shocks in a weak spiral potential (Wada & Koda 2004;Khoperskov et al. 2011). Once we select stars in the guiding space features (see Fig. A.1g) we calculate the radial velocity distribution (see Fig. A.1h). For this model and for the observed SDS structures, error bars are determined using the standard deviation of the measurements. The distributions show that the stars associated with spiral arms move inwards (negative radial velocity). This systematic motion explains the mean radial velocity of stars in the spiral arms and their skewed distribution profile (see Fig. A.1h). This result is further evidence that stellar density structures with similar kinematics (SDS 1,3,4,6) correspond to the spiral arms of the Milky Way. In Model I, the spiral structure is a slowly rotating pattern with a corotation outside 10 kpc, where a systematic outward motion (positive radial velocity) of stars is observed in agreement with the density wave theory.

A.2. Model II: barred galaxy model
Model II is a purely collisionless N-body simulation of a composite disc galaxy with live dark matter halo with a total number of particles of 10 7 for stars and 0.5 × 10 7 for dark matter. For the N-body system integration, we used our parallel version of the TREE-GRAPE code with multi-thread usage under the SSE and AVX instructions. The tolerance parameter of the tree-code is Θ = 0.7, and for smoothing we use a Plummer potential with 50 pc. For the time integration, we used a leapfrog integrator with a fixed step size of 0.2 Myr. This simulation was previously used to reproduce the morphology of the metal-rich and metal-poor stellar populations in the Milky Way bulge (Fragkoudi et al. 2017(Fragkoudi et al. , 2018, to understand the origin of the phase-space spirals in isolated Milky Way-type galaxies (Khoperskov et al. 2019), and to explain the ridges in the R−V φ plane of the Milky Way (Martinez-Medina et al. 2019;Fragkoudi et al. 2019). In this work we refer in our analysis to the snapshot at 4 Gyr from the beginning of evolution, when no prominent spiral structure is seen, and the bar is the dominant non-axisymmetric structure perturbing the phasespace of the simulated galaxy (see Fig. A.2a).
In Fig. A.2 we show the density distribution of stars in both spatial and guiding space where two prominent quasi-circular distinct features appear in the guiding space. The motion of a star in a galactic plane is described by radial (κ -epicyclic frequency) and angular oscillation (Ω is the rotation around the galactic centre). To confirm that these overdensities are related to the bar resonances we measure the orbital frequencies (Ω, κ) (Binney & Spergel 1982;Ceverino & Klypin 2007;Halle et al. 2018) of all particles over a 1 Gyr period around 4 Gyr. Each of the bar resonances is characterized by a certain value of F = (Ω − Ω bar )/κ (e.g. for the bar corotation F = 0, and for inner (ILR) and outer Lindblad (OLR) resonances F = −0.5 and 0.5, respectively). In Fig. A.2c the overall distribution peaks at the ILR, which is expected for a galaxy with a strong bar, while the OLR is rather weak but slightly stronger than the corotation. The coloured lines show the frequency ratios for stars associated with two prominent features in the guiding phase-space where the inner ring-like structure frequency ratio peaks at corotation (blue) and the outer ring consists in the OLR stars (red). Therefore, the two features we found in the guiding space of the barred galaxy correspond to the corotation and the OLR, or, in other words, the main bar resonances are imprinted in the guiding space as large-scale overdensities. Stellar density distributions and radial velocity maps of stars at the main bar resonances in guiding and spatial scales are presented in Fig. A.2d-h L8, page 6 of 9 where stars that are trapped at the bar resonances cover a large area across the disc. The radial velocity distribution along the bar corotation region is roughly symmetric, and the variations in the mean value are caused by sequential gain and loss of angular momentum, while there are no net changes in angular momentum over a long period of time (Ceverino & Klypin 2007). At the bar OLR stars librate around the 2 : 1 closed periodic orbits that correspond to the x 1 (1) and x 1 (2) families, which contribute to creating the asymmetry in radial velocity along the OLR (Fragkoudi et al. 2019). These changes in the radial velocities are imprinted in the radial velocity distribution for stars at the resonances (see Fig. A.2i); in particular, for a disc region similar to that available in the Milky Way (|X| < 5 kpc and Y < 0), the mean radial velocity is negative and the distribution spreads over a wide range of values. At the OLR, the radial velocity distribution is dominated by stars on x 1 (1) orbits L8, page 7 of 9 A&A 634, L8 (2020) density in the guiding space. Two distinct quasi-circular overdensities appear close to the end of the bar and in the outer disc indicated by blue and red circles, respectively. Panel c: spectral analysis of orbits of stars where the fraction of stars with different F = (Ω − Ω bar )/κ ratio is shown for all disc stars (black) and stars between the blue and red circles in frame (b). Frequency ratios of −0.5, 0, and 0.5 correspond to the main resonances of the bar, the outer Lindblad, corotation, and inner Lindblad resonances, respectively. Panel d: guiding space density distribution of stars at the bar corotation (|F| < 0.05) and at the OLR (|F + 0.5| < 0.05). Panel e: density distribution of stars corotating with the bar in spatial coordinates. Panel f: density distribution of stars at the bar OLR in spatial coordinates. Panel g: mean radial velocity distribution of all stars in the guiding space. Panel h: mean radial velocity distribution of stars at the bar corotation and the OLR in the guiding space. Panel i: distributions of the radial velocities of stars at the bar corotation (blue) and at the OLR (red).
which have positive radial velocities and, on average, this OLR patch moves outwards (positive radial velocity). In Fig. A.2i the error bars in the radial velocity distributions are from comparing 20 simulation snapshots within a ±0.5 Gyr range, using the standard deviation of the measurements. The phase-space and kinematical features that we found for stars at the resonances are in qualitative agreement with those we discovered in the Milky Way disc by analysing the Gaia DR2 (see Fig. 2). Therefore, we confirm that the SDS 2 is likely to be associated with the corotation, while the SDS 5 corresponds to the OLR of the bar.