# A dynamically young and perturbed Milky Way disk

T. Antoja, A. Helmi, M. Romero-Gomez, D. Katz, C. Babusiaux, R. Drimmel, D. Evans, F. Figueras, E. Poggio, C. Reylé, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

T. Antoja, A. Helmi, M. Romero-Gomez, D. Katz, C. Babusiaux, et al., A dynamically young and perturbed Milky Way disk. Nature, 2018, 561 (7723), pp.360 - 362. 10.1038/s41586-018-0510-7. hal-01878376

HAL Id: hal-01878376

https://hal.science/hal-01878376v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A dynamically young and perturbed Milky Way disk

T. Antojá<fnref n="1" />*, A. Helmi², M. Romero-Gómez<fnref n="1" />, D. Katz³, C. Babusiaux³,⁴, R. Drimmel⁵, D. W. Evans⁶, F. Figueras<fnref n="1" />, E. Poggio⁵,⁷, C. Reylé⁸, A. C. Robin⁸, G. Seabroke⁹ & C. Soubiran¹⁰

The evolution of the Milky Way disk, which contains most of the stars in the Galaxy, is affected by several phenomena. For example, the bar and the spiral arms of the Milky Way induce radial migration of stars<fnref n="1" /> and can trap or scatter stars close to orbital resonances². External perturbations from satellite galaxies can also have a role, causing dynamical heating of the Galaxy³, ring-like structures in the disk⁴ and correlations between different components of the stellar velocity⁵. These perturbations can also cause ‘phase wrapping’ signatures in the disk⁶⁻⁹, such as arched velocity structures in the motions of stars in the Galactic plane. Some manifestations of these dynamical processes have already been detected, including kinematic substructure in samples of nearby stars¹⁰⁻¹², density asymmetries and velocities across the Galactic disk that differ from the axisymmetric and equilibrium expectations¹³, especially in the vertical direction¹¹,¹⁴⁻¹⁶, and signatures of incomplete phase mixing in the disk⁷,¹²,¹⁷,¹⁸. Here we report an analysis of the motions of six million stars in the Milky Way disk. We show that the phase-space distribution contains different substructures with various morphologies, such as snail shells and ridges, when spatial and velocity coordinates are combined. We infer that the disk must have been perturbed between 300 million and 900 million years ago, consistent with estimates of the previous pericentric passage of the Sagittarius dwarf galaxy. Our findings show that the Galactic disk is dynamically young and that modelling it as time-independent and axisymmetric is incorrect.

Gaia is a European Space Agency (ESA) mission that was designed primarily to investigate the origin, evolution and structure of the Milky Way, and has recently delivered the largest and most precise census of positions, velocities and other stellar properties of more than a billion stars. By exploring the phase space (positions and velocities) of more than six million stars within a few kiloparsecs of the Sun in the Galactic disk from Gaia data release 2 (DR2; see Methods)¹⁹, we find that certain

phase-space projections (Figs. 1a, 2) have many substructures that had not been predicted by existing models. These substructures had remained blurred until now, owing to the limitations on the number of stars in and the precision of previously available datasets.

In Fig. 1a we show the projection of phase space in vertical position and velocity, $Z-V_Z$. The stars follow a curled, spiral-shaped distribution, the density of which increases towards the leading edge of the spiral. Figure 1b, c demonstrates that this ‘snail shell’ pattern is still present when the stars are colour-coded according to their radial and azimuthal velocities, $V_R$ and $V_\phi$, which implies a strong correlation between the vertical and in-plane motions of the stars. The pattern is particularly pronounced in the case of $V_\phi$ (Fig. 1c), even up to $V_Z \approx 40 \text{ km s}^{-1}$. Furthermore, we see a gradient of different values of $V_\phi$ across the spiral shape, which follows the density variations. Details of the relationship between the snail shell and the other velocity features observed in the solar neighbourhood are shown in Extended Data Fig. 1.

The spiral shape in Fig. 1a is clearly reminiscent of the effects of phase mixing in two dimensions that have been discussed in several areas of astrophysics²⁰⁻²² and in quantum physics²³, but never in the context of dynamical models of the Galactic disk. This process can be better understood by using a toy model. Consider a Galaxy model in which the vertical potential of the Galaxy can be approximated by an anharmonic oscillator (see equation (1) in Methods). In this approximation, the vertical frequencies of oscillation $\nu$ depend on the amplitude of the oscillation A and the Galactocentric radius R, to first order²² (see equation (2) in Methods). By assuming that stars follow a simple harmonic oscillation with these frequencies, their movement over time $t$ is described by $Z = A\cos[\phi(t)]$ and $V_Z = -A\nu\sin[\phi(t)]$ (where $\phi(t)$ is the orbital phase), which traces an oval shape in the clockwise direction in the $Z-V_Z$ projection. However, stars revolve at different angular speeds depending on their frequency. Thus, an ensemble of stars will stretch out in phase space, with the range of frequencies causing a spiral

<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Vertical positions and velocities of the stars. The plots show the distribution of stars in the vertical position–velocity ($Z-V_Z$) plane from our sample of Gaia data for stars with Galactocentric radii of 8.24 kpc $< R < 8.44$ kpc. **a**, Two-dimensional histogram in bins of $\Delta Z = 0.02$ kpc and $\Delta V_Z = 1$ km s$^{-1}$, with the darkness of the colour scale</figcaption></figure>

proportional to the number of stars. **b**, $Z-V_Z$ plane coloured as a function of median radial velocity $V_R$ in bins of $\Delta Z = 0.02$ kpc and $\Delta V_Z = 1$ km s$^{-1}$. **c**, Same as **b**, but for the azimuthal velocity $V_\phi$. $V_R$ and $V_\phi$ are positive towards the Galactic anticentre and the direction of Galactic rotation, respectively.

<footnote n="1">¹Institut de Ciències del Cosmos, Universitat de Barcelona (IEEC-UB), Barcelona, Spain. ²Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands. ³GEPI, Observatoire de Paris, Université PSL, CNRS, Meudon, France. ⁴Université Grenoble Alpes, CNRS, IPAG, Grenoble, France. ⁵INAF—Osservatorio Astrofisico di Torino, Pino Torinese, Italy. ⁶Institute of Astronomy, University of Cambridge, Cambridge, UK. ⁷Université di Torino, Dipartimento di Fisica, Torino, Italy. ⁸Institut UTINAM, CNRS UMR6213, Université Bourgogne Franche-Comté, OSU THETA Franche-Comté Bourgogne, Observatoire de Besançon, Besançon, France. ⁹Mullard Space Science Laboratory, University College London, Dorking, UK. ¹⁰Laboratoire d'astrophysique de Bordeaux, Université Bordeaux, CNRS, Pessac, France. *e-mail: tantoja@fqa.ub.edu</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2 | Positions and velocities of the stars in the disk plane.**
Distribution of azimuthal velocities $V_\phi$ as a function of Galactocentric radius $R$ for all stars in our sample with six-dimensional phase space-coordinates from Gaia DR2, shown as a two-dimensional histogram in bins of $\Delta V_\phi = 1 \text{ km s}^{-1}$ and $\Delta R = 0.01 \text{ kpc}$.</figcaption></figure>

shape in this projection. The detailed time evolution of stars in this toy model is described in Methods and shown in Extended Data Fig. 3. As time goes by, the spiral gets more tightly wound and eventually this process of phase mixing leads to a spiral that is so tightly wound that the coarse-grained distribution appears to be smooth. The clarity of the spiral shape in the $Z-V_Z$ plane revealed by Gaia DR2 implies that this time has not yet been reached in the Milky Way, providing evidence that phase mixing is currently taking place in the Galactic disk.

This interpretation also implies that the shape of the spiral can be used to obtain information about: (i) the shape of the potential, which determines the vertical frequencies; (ii) the starting time of phase mixing; and (iii) the type of perturbation that brought the disk into a non-equilibrium state, which sets the initial conditions for the phase-mixing event that we are witnessing. For instance, we can estimate the time $t$ of the event from the separation between two consecutive spiral turns, because these have a phase separation of $(\nu_2 t + \phi_0) - (\nu_1 t + \phi_0) = 2\pi$, where the subscripts ‘1’ and ‘2’ indicate two consecutive turns of the spiral. Therefore, assuming that the initial phase $\phi_0$ is the same for turns 1 and 2, we have $t = 2\pi/(\nu_2 - \nu_1)$. Using several potentials for the Milky Way (see Methods), we estimate that the vertical phase-mixing event started between 300 Myr and 900 Myr ago. The toy model that illustrates this process is shown in Fig. 3a, which depicts a snail shell that formed after 500 Myr from an ensemble of stars with a starting distribution that is out of equilibrium; this modelled snail shell is similar to the one seen in the data.

A possible perturbation that might have initiated the on-going vertical phase mixing that we observe is the influence of a satellite galaxy. In particular, the last pericentre of the orbit of the Sagittarius dwarf galaxy has been shown to have strong effects on the stellar disk$^{4,8,9}$. In addition, most models place this pericentric passage between 200 Myr and 1,000 Myr ago$^{9,24,25}$, which is consistent with our findings. Nevertheless, other processes that may induce snail shell patterns include the formation of the central bar and of the transient spiral structure, provided that these processes can induce vertical asymmetries, other global changes in the potential, and the dissolution of a massive stellar system such as a cluster or accreted satellite.

Another phase-space projection in which the Gaia data have a markedly different appearance is that in azimuthal velocity and cylindrical radius, $V_\phi - R$ (Fig. 2). Although this phase-space projection has been explored previously with other data$^{26}$, the spatial coverage, high sampling and unprecedented precision of the Gaia data reveal many thin diagonal ridges that were not previously evident. The arches in the velocity-space projection $V_R - V_\phi$ in the solar neighbourhood that were discovered recently in Gaia data$^{12}$ (see Extended Data Fig. 1a) are projections of these diagonal ridges, but at a fixed Galactic position. Therefore, Fig. 2 reveals that arched velocity structures must be present at many different radii, but have thus far been unexplored, and that their characteristics vary with distance from the Galactic centre,

<figure><img src="image_5.png" /><figcaption>**Fig. 3 | Models of the phase-space distribution of the Galaxy disk.**
a, Modelled spiral shape created in the vertical position–velocity ($Z-V_Z$) plane as a result of phase mixing in the evolution of an ensemble of particles for 500 Myr in a Galactic potential, starting from a distribution that is out of equilibrium, presumably after a perturbation. b, Modelled diagonal ridges created in the distribution of azimuthal velocities $V_\phi$ as a function of Galactocentric radius $R$ as a result of the phase mixing in the evolution of an ensemble of particles for 1,000 Myr in a Galactic potential, starting from a distribution that is out of equilibrium. c, Same as b, but for diagonal ridges created as a result of the effects of the barred potential and its resonant structure. See Methods for more details.</figcaption></figure>

diminishing their velocity towards the outskirts of the Galaxy in a continuous way.

These diagonal ridges could be signatures of phase mixing in the horizontal direction, as has been predicted for the arches in velocity space$^{6–8}$. Alternatively, the bar and the spiral arms could induce diagonal ridges through their resonant orbital structure, creating regions in phase space of stable and unstable orbits$^{27}$, and hence with over-densities and gaps. The toy model of phase mixing (Fig. 3b) and a disk simulation with a Galactic potential that contains a bar (Fig. 3c) both show several diagonal ridges. The $V_\phi$ separation of consecutive ridges in the data is about 10 km s$^{-1}$. Comparing this separation to that of our toy model (see Methods) indicates that, if these ridges are caused by phase mixing from a single perturbation, it should have taken place longer ago than the perturbation that gave rise to the vertical mixing. This is consistent with the timing derived using the separation between arches in the local velocity plane$^7$ and with the existence of a group of co-moving stars that appear not to be fully phase mixed vertically$^{18}$, which suggests that another perturbation occurred about 2 Gyr ago. The relationship between the various features is not clear, and it
is not unlikely that there are or were several perturbations creating
superposed features.

Our interpretation of the features that we found is based on toy
models, the main limitations of which are their lack of self-consistency,
the choice of initial conditions not necessarily reflecting those stem-
ming from the impact of a satellite galaxy, and the fact that we study
separately both the effects of resonances and phase mixing, and the
different phase-space dimensions (horizontal and vertical) involved.
A challenging task for the future will be to model our findings taking
into account collective effects, such as in the perturbative regime²⁸, and
with self-consistent N-body models²⁵.

Online content

Any methods, additional references, Nature Research reporting summaries, source data, statements of data availability and associated accession codes are available at
https://doi.org/10.1038/s41586-018-0510-7.

Received: 26 April; Accepted: 26 July 2018;
Published online 19 September 2018.

1. Sellwood, J. A. & Binney, J. J. Radial mixing in galactic discs. *Mon. Not. R. Astron. Soc.* **336**, 785–796 (2002).

2. Contopoulos, G. & Grosbol, P. Stellar dynamics of spiral galaxies: nonlinear effects at the 4/1 resonance. *Astron. Astrophys.* **155**, 11–23 (1986).

3. Quinn, P. J., Hernquist, L. & Fullagar, D. P. Heating of galactic disks by mergers. *Astrophys. J.* **403**, 74–93 (1993).

4. Purcell, C. W., Bullock, J. S., Tollerud, E. J., Rocha, M. & Chakrabarti, S. The Sagittarius impact as an architect of spirality and outer rings in the Milky Way. *Nature* **477**, 301–303 (2011).

5. D'Onghia, E., Madau, P., Vera-Ciro, C., Quillen, A. & Hernquist, L. Excitation of coupled stellar motions in the Galactic disk by orbiting satellites. *Astrophys. J.* **823**, 4 (2016).

6. Fux, R. Order and chaos in the local disc stellar kinematics induced by the Galactic bar. *Astron. Astrophys.* **373**, 511–535 (2001).

7. Minchev, I. et al. Is the Milky Way ringing? The hunt for high-velocity streams. *Mon. Not. R. Astron. Soc.* **396**, L56–L60 (2009).

8. Gómez, F. A., Minchev, I., Villalobos, Á., O'Shea, B. W. & Williams, M. E. K. Signatures of minor mergers in Milky Way like disc kinematics: ringing revisited. *Mon. Not. R. Astron. Soc.* **419**, 2163–2172 (2012).

9. de la Vega, A., Quillen, A. C., Carlin, J. L., Chakrabarti, S. & D'Onghia, E. Phase wrapping of epicyclic perturbations in the Wobbly galaxy. *Mon. Not. R. Astron. Soc.* **454**, 933–945 (2015).

10. Eggen, O. J. Star streams and Galactic structure. *Astron. J.* **112**, 1595–1613 (1996).

11. Dehnen, W. The distribution of nearby stars in velocity space inferred from HIPPARCOS data. *Astron. J.* **115**, 2384–2396 (1998).

12. Gaia Collaboration. Gaia data release 2: mapping the Milky Way disc kinematics. *Astron. Astrophys.* **616**, A11 (2018).

13. Siebert, A. et al. Detection of a radial velocity gradient in the extended local disc with RAVE. *Mon. Not. R. Astron. Soc.* **412**, 2026–2032 (2011).

14. Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S. & Chen, H.-Y. Galactoseismology: discovery of vertical waves in the Galactic disk. *Astrophys. J.* **750**, L41 (2012).

15. Schönrich, R. & Dehnen, W. Warp, waves, and wrinkles in the Milky Way. *Mon. Not. R. Astron. Soc.* **478**, 3809–3824 (2018).

16. Quillen, A. C. et al. The GALAH survey: stellar streams and how stellar velocity distributions vary with Galactic longitude, hemisphere and metallicity. *Mon. Not. R. Astron. Soc.* **478**, 228–254 (2018).

17. Gómez, F. A. et al. Signatures of minor mergers in the Milky Way disc – I. The SEGUE stellar sample. *Mon. Not. R. Astron. Soc.* **423**, 3727–3739 (2012).

18. Monari, G., et al. Coma Berenices: the first evidence for incomplete vertical phase-mixing in local velocity space with RAVE—confirmed with Gaia DR2. *Res. Notes AAS* **2**, 32 (2018).

19. Gaia Collaboration. Gaia data release 2: summary of the contents and survey properties. *Astron. Astrophys.* **616**, A1 (2018).

20. Tremaine, S. The geometry of phase mixing. *Mon. Not. R. Astron. Soc.* **307**, 877–883 (1999).

21. Afshordi, N., Mohayaee, R. & Bertschinger, E. Hierarchical phase space structure of dark matter haloes: tidal debris, caustics, and dark matter annihilation. *Phys. Rev. D* **79**, 083526 (2009).

22. Candlish, G. N. et al. Phase mixing due to the Galactic potential: steps in the position and velocity distributions of popped star clusters. *Mon. Not. R. Astron. Soc.* **437**, 3702–3717 (2014).

23. Manfredi, G. & Feix, R. M. Theory and simulation of classical and quantum echoes. *Phys. Rev. E* **53**, 6460–6470 (1996).

24. Law, D. R. & Majewski, S. R. The Sagittarius dwarf galaxy: a model for evolution in a triaxial Milky Way halo. *Astrophys. J.* **714**, 229–254 (2010).

25. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N. & Besla, G. The influence of Sagittarius and the Large Magellanic Cloud on the Milky Way galaxy. *Mon. Not. R. Astron. Soc.* https://doi.org/10.1093/mnras/sty1574 (2018).

26. Monari, G., Kawata, D., Hunt, J. A. S. & Famaey, B. Tracing the Hercules stream with Gaia and LAMOST: new evidence for a fast bar in the Milky Way. *Mon. Not. R. Astron. Soc.* **466**, L113–L117 (2017).

27. Michtchenko, T. A., Lépine, J. R. D., Barros, D. A. & Vieira, R. S. S. Combined dynamical effects of the bar and spiral arms in a Galaxy model. Application to the solar neighbourhood. *Astron. Astrophys.* **615**, A10 (2018).

28. Fouvry, J.-B., Binney, J. & Pichon, C. Self-gravity, resonances, and orbital diffusion in stellar disks. *Astrophys. J.* **806**, 117 (2015).

**Acknowledgements** This work made use of data from ESA mission Gaia (https://www.cosmos.esa.int/gaia), which was processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This project received funding from the European Union's Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement number 745617. This work was supported by the MDM-2014-0369 of ICCUB (Unidad de Excelencia 'María de Maeztu') and the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement GENIUS FP7-606740. A.H. acknowledges financial support from a VICI grant from the Netherlands Organisation for Scientific Research (NWO). We acknowledge the MINECO (Spanish Ministry of Economy) through grants ESP2016-80079-C2-1-R (MINECO/FEDER, UE) and ESP2014-55996-C2-1-R (MINECO/FEDER, UE). This work been funded in part by the Agenzia Spaziale Italiana (ASI) through contract 2014-025-R.1.2015 through the Italian Istituto Nazionale di Astrofisica (INAF). E.P. acknowledges the financial support of the 2014 PhD fellowship programme of INAF.

**Author contributions** T.A. contributed to the sample preparation, analysed and interpreted the data, performed most of the modelling and wrote the paper together with A.H. A.H. also provided interpretation of the findings. M.R.-G. performed the simulation with the barred potential and contributed to sample preparation. D.K., C.B., R.D., D.W.E., F.F., E.P., C.R., A.C.R., G.S. and C.S. contributed to sample preparation and validation of the Gaia data. All authors reviewed the manuscript.

**Competing interests**: The authors declare no competing interests.

**Additional information**

*Extended data* is available for this paper at https://doi.org/10.1038/s41586-018-0510-7.

**Reprints and permissions information** is available at http://www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to T.A.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
## METHODS

**Data and samples selection.** We used Gaia DR2 sources for which the six-dimensional phase-space coordinates can be computed, that is all stars for which a five-parameters astrometric solution (sky positions, parallax and proper motions) and a line-of-sight velocity are available. We selected only stars with positive parallaxes $\bar{\omega}$ with relative uncertainty less than 20%, that is, satisfying $\bar{\omega}/\sigma_{\bar{\omega}} > 5$. This selection ensures that $1/\bar{\omega}$ is a reasonably good estimator of the distance to the stars$^{29}$; alternatively, we also used Bayesian distances (see below). This sample has 6,376,803 stars and has been well studied and characterized previously$^{12}$. The data were obtained directly through the following query in the public Gaia archive (https://gea.esac.esa.int/archive):

```sql
SELECT G.source_id, G.radialVelocity, G.radialVelocity_error, G.r_a, G.r_a_error, G.dec, G.dec_error, G.parallax, G.parallax_error, G.pmra, G.pmra_error, G.pmdec, G.pmdec_error, G.r_a_dec_corr, G.r_a_parallax_corr, G.r_apmdec_corr, G.r_apmdec_corr, G.dec_parallax_corr, G.decpmra_corr, G.dec_pmddec_corr, G.parallaxpmra_corr, G.parallax_pmddec_corr, G.pmra_pmddec_corr
```

```sql
FROM gaiadr2.gaia_source G
```

```sql
WHERE G.radialVelocity IS NOT Null AND G.parallax_over_error>5.
```

From the five-parameter astrometric solution and line-of-sight velocities ($\alpha$, $\delta$, $\bar{\omega}$, $\mu_{\alpha}^*$, $\mu_{\delta}$, $V_{10s}$) of these stars, we derived distances (as $1/\bar{\omega}$), positions and velocities in the cylindrical Galactic reference frame ($R$, $\phi$, $Z$, $V_R$, $V_{\phi}$, $V_Z$). For convenience, we took $\phi$ to be positive in the direction of Galactic rotation and with the origin at the Sun-Galactic centre line. For these transformations, we adopted a vertical distance of the Sun above the plane of$^{30}$ 27 pc, a distance of the Sun to the Galactic centre$^{31}$ of $R_{\odot} = 8.34$ kpc and a circular velocity at the Sun radius of$^{31}$ $V_C(R_{\odot}) = 240$ km s$^{-1}$. We assumed a peculiar velocity of the Sun with respect to the local standard of rest of$^{32}$ $(U_{\odot}, V_{\odot}, W_{\odot}) = (11.1, 12.24, 7.25)$ km s$^{-1}$. Our choice of values gives $(V_C(R_{\odot}) + V_{\odot})/R_{\odot} = 30.2$ km s$^{-1}$ kpc$^{-1}$, which is compatible with the reflex motion of Sgr A*$^{33}$. To derive the uncertainties in these coordinates, we propagate the full covariance matrix. The median uncertainties in $V_R$, $V_{\phi}$ and $V_Z$ are 1.4 km s$^{-1}$, 1.5 km s$^{-1}$ and 1.0 km s$^{-1}$, respectively, and 80% of stars have uncertainty smaller than 3.3 km s$^{-1}$, 3.7 km s$^{-1}$ and 2.2 km s$^{-1}$, respectively, in these velocities. The positions in the Cartesian coordinates X-Y and X-Z of the sample are shown in Extended Data Fig. 2.

For part of our study, we selected from our sample the 935,590 stars located in the solar Galactic cylindrical ring, with Galactocentric radius 8.24 kpc $< R < 8.44$ kpc (dotted lines in Extended Data Fig. 2). For this selection, the median uncertainties in $V_R$, $V_{\phi}$ and $V_Z$ are 0.5 km s$^{-1}$, 0.8 km s$^{-1}$ and 0.6 km s$^{-1}$, respectively, and 80% of stars have uncertainty smaller than 1.1 km s$^{-1}$, 2.0 km s$^{-1}$ and 1.3 km s$^{-1}$, respectively, in these velocities.

The velocity uncertainties are significantly smaller than the sizes of the substructures detected; together with the large number of stars in our samples, this is what made their detection possible. Although there are some correlations between the astrometric Gaia observables$^{34}$, these are not responsible for the correlations and substructures seen in our phase-space plots. This is because the stars in our sample are distributed across all sky directions, and the phase-space coordinates come from combinations of astrometric measurements and radial velocities in different contributions depending on the direction on the sky. Besides, the astrometric correlations for our sample are small (less than 0.2 in their absolute value for more than 50% of stars); this, combined with the small errors, makes their contribution insignificant.

Alternatively, we used distances determined through a Bayesian inference method using the existing implementation in TOPCAT$^{35}$, taking the mode of the posterior distribution and a prior of an exponentially decreasing density of stars with scale length of$^{36}$ 1.35 kpc. We found that the differences between this distance determination and the inverse of the parallax are between -2% and 0.6% for 90% of the 6,376,803 stars with $\bar{\omega}/\sigma_{\bar{\omega}} > 5$, which is expected for small relative errors in parallax. Consequently, the phase-space diagrams presented here vary only at the pixel level. The diagrams do not vary even when using the set of 7,183,262 of stars with available radial velocities, which includes stars with larger parallax errors and with negative parallaxes, for which the estimator of the inverse of the parallax would yield unphysical distances. When using another alternative set of Bayesian distances specifically derived for stars from Gaia DR2 with radial velocities using a different prior$^{37}$, we found the differences between these distances and the inverse of the parallax to be between -9% and 5% for 90% of the stars—slightly larger than before, but again with no noticeable effects on the phase-space projections examined here.

**Models for vertical phase mixing.** We first reproduced the spiral shape observed in the Z-VZ plane with Gaia data by using a toy model. The classic harmonic oscillator is often used to describe the vertical movement of stars in galactic disks under the epicyclic theory$^{38}$. However, in this approximation, which is valid only

for very small-amplitude orbits for which the potential changes little vertically, stars have the same vertical oscillatory frequency $\nu$ and there is no phase mixing, unless orbits at different guiding radius, and therefore with different frequencies, are considered. Instead, we use an anharmonic oscillator with potential

$$ \Phi(Z) \propto -\alpha_0 + \frac{1}{2}\alpha_1 Z^2 - \frac{1}{4}\alpha_2 Z^4 \quad (1) $$

We adopt coefficients $\alpha_0$, $\alpha_1$ and $\alpha_2$ that correspond to the expansion for small Z (derived elsewhere$^{22}$) of a Miyamoto-Nagai potential$^{39}$ with $a=6.5$ kpc, $b=0.26$ kpc and $M=10^{11} M_{\odot}$. These coefficients depend on Galactocentric radius $R$ because the vertical pull depends on the distance to the Galactic centre. In this anharmonic potential, the frequencies of oscillation are described by

$$ \nu(A, R) = \alpha_1(R)^{1/2} \left[ 1 - \frac{3\alpha_2(R)A^2}{8\alpha_1(R)} \right] \quad (2) $$

where $\nu_0 = \alpha_1^{1/2}$ is the vertical frequency in the epicyclic approximation.

Given an initial distribution of stars with $Z(t=0)$ and $V_Z(t=0)$, the vertical amplitudes of the orbits can be derived through the conservation of energy and using the fact that at the vertical turn-around point of the orbit ($V_Z=0$) the (vertical) kinetic energy is null$^{22}$. Assuming that stars follow a simple harmonic oscillation (but with different frequencies), the movement of the stars with time is described by

$$ \begin{aligned} Z &= A \cos[\nu(A, R)t + \phi_0] \\ V_Z &= -A\nu(A, R)\sin[\nu(A, R)t + \phi_0] \end{aligned} \quad (3) $$

where the initial phase of the stars $\phi_0 = \phi(t=0)$ is obtained from the initial distribution of Z and $V_Z$ and the corresponding amplitudes.

The phase-space evolution described above is shown in Extended Data Fig. 3a-c. Initially, the particles follow a Gaussian distribution in $Z(t=0)$ and $V_Z(t=0)$ with mean and dispersion of -0.1 kpc and 0.04 kpc, and -2 km s$^{-1}$ and 1 km s$^{-1}$, respectively. We located all particles at the same Galactocentric radius of $R = 8.5$ kpc; therefore, they all move under the same functional form of the vertical potential. The initial conditions are shown in Extended Data Fig. 3a, in which we have colour-coded the particles according to their period. Following equation (3), each star follows a clockwise rotation in the Z-VZ plane. However, they do so at different angular speeds: stars with smaller period located at closer distances from the mid-plane ($Z=0$) revolve faster than those located at largest distances from the mid-plane. The range of frequencies is therefore what creates the spiral shape. Extended Data Fig. 3b shows the evolution of the system for three initial phases of the time evolution when the spiral shape begins to form. Extended Data Fig. 3c shows the spiral shape after 1,000 Myr of evolution.

In the Gaia data (Fig. 1), we do not see a thin spiral but a thick one, with many of the stars in the volume participating in it. A similar effect was observed with our toy model when we included particles at different radius for which the vertical potential and the range of amplitudes and frequencies changes. In Extended Data Fig. 3d-f we let a similar system evolve as in Extended Data Fig. 3a-c but starting with the initial radius following a skewed normal distribution, which creates a density that decreases with radius, as in galaxy disks, with a skewness of 10, location parameter of 8.4 kpc and scale parameter of 0.2 kpc. The spiral structure is now thickened, similarly to the data, with a higher density of stars at the leading edge of the spiral.

To estimate the time of the phase-mixing event from the spiral seen in the Gaia data (Fig. 1) using

$$ t = \frac{2\pi}{\nu_2 - \nu_1} \quad (4) $$

we need to locate two consecutive turns of the spiral and estimate their vertical frequencies from their amplitudes and mean radius. For this, we used Extended Data Fig. 5, which is colour-coded as a function of median guiding radius. The guiding radius is approximated as $R_g \approx V_{\phi}R_{\odot}/[V_C(R_{\odot})]$, under the hypothesis of a flat rotation curve, where we used the values of $R_{\odot} = 8.34$ kpc and $V_C(R_{\odot}) = 240$ km s$^{-1}$ assumed in the coordinate transformation of the data. From Extended Data Fig. 5 we see that the density gradient across the spiral shape is created by stars with different guiding radius that arrive at the solar neighbourhood owing to their different amplitudes of (horizontal) radial oscillation. To determine two consecutive turns of the spiral, we focused on stars at the turn-around points ($V_Z=0$) near the leading edges of the spiral. By visual inspection, we determined an approximate range of Z in which the turn-around points are located in Extended Data Fig. 5, concentrating on red colours, for which the spiral is well defined. The ranges of the turn-around points are marked with vertical lines and listed, together with the middle value, in Extended Data Table 1. For these turn-around points, the
amplitudes are simply $A = Z$ and from the colour bar we note that the median $R_g$ is around 8.2 kpc. Small changes in this value do not change substantially our final determination of the time of the perturbation.

To estimate the vertical orbital frequencies of these turn-around points, we could not use the toy model presented above because it is valid only for oscillations with small amplitude A—in particular, smaller than the vertical scale b of the potential ($A \ll 0.26$ kpc). Therefore, we took an existing model⁴⁰ with updated parameters that fit current estimates such as for the Sun Galactocentric radius and the circular velocity curve⁴¹. We computed the vertical frequency numerically in a grid of different radii and vertical amplitudes by integrating orbits and measuring their vertical periods (Extended Data Fig. 4a). The vertical frequency can change along the orbits for stars with large eccentricities in the horizontal direction, but here for simplicity we put all particles on nearly circular orbits. We estimated the vertical frequency at each turn-around position directly by interpolating the calculations of Extended Data Fig. 4a using the estimated values for the amplitude and radius.

Finally, taking each pair of turning points, we obtained an estimate of the time since the perturbation using equation (4). As an example, the two turning points (amplitudes) of the left part of the spiral are located at $-0.59 \pm 5$ kpc and $-0.23 \pm 5$ kpc, respectively. These correspond to vertical frequencies of $0.058 \pm 0.002$ rad Myr$^{-1}$ and $0.072 \pm 0.002$ rad Myr$^{-1}$ for $R_g = 8.2$ kpc, which gives a time of $461^{+183}_{-105}$ Myr. We repeated the same procedure for the second pair of consecutive turning points, and also for the edges of the spiral at $Z=0$ (mid-plane points, which have $V_Z = A/r$), estimating the frequencies by interpolating in Extended Data Fig. 4b. The mid-plane positions are marked as horizontal lines in Extended Data Fig. 5. All results are summarized in Extended Data Tables 1 and 2. The mean of the three time estimates is 510 Myr and the minimum and maximum times from the uncertainty ranges are 356 Myr and 856 Myr.

We tested the dependence of our time determination on the potential model used by using a different model⁴². Compared to our previous model, this one has different shapes for the halo, disks and bulge, a different total mass, and includes thin and thick disks and two gas disks. The frequencies in this model are smaller by 4% on average and smaller than 6% for 90% of the points in the grid of Extended Data Fig. 4a. By repeating the whole process to determine the perturbation time, we obtained 528 Myr with minimum and maximum times from the uncertainty ranges of 361 Myr and 899 Myr, very similar to our previous determination.

Our determination is subject to several approximations: (1) we used the vertical frequencies of orbits with conditions of circularity on the plane and for certain assumed Galactic potential models; (2) we considered a unique guiding radius; and (3) we took equal initial phases for the turn-around and mid-plane points.

We finally ran a simulation (Fig. 3a) by integrating 100,000 test-particle orbits in the updated model⁴¹ with initial vertical positions and velocities following Gaussian distributions centred at $Z = -0.4$ kpc and $V_Z = -5$ km s$^{-1}$ and with dispersions of 0.15 kpc and 2 km s$^{-1}$, respectively. Horizontally in the disk plane, the test particles were distributed following a skewed normal distribution in radius $R$, with a scale parameter of 0.8 kpc, location parameter of 8, skewness of 10, and all particles at an azimuthal angle $\phi = 0$. The horizontal velocities were set to 0 for the radial component and to the circular velocity at the particle's radius and $Z=0$ for the azimuthal component. These are initial conditions of circularity on the Galactic plane but not necessarily for orbits with large excursions in Z. The particle's orbits were integrated forwards in time for 500 Myr as estimated from the data. This is not meant to be a fit to the data because we have not explored all possible initial configurations that could lead to a similar spiral shape. However, we see that an initial distribution that is asymmetric in Z, with most particles located at positive or negative Z, is required to obtain a single spiral instead of a symmetrical double one. **Model for the horizontal phase mixing.** We used a toy model to reproduce the diagonal ridges observed in the $V_\phi-R$ plane (Fig. 3b). This model is built by integrating orbits in the Galactic potential of ref.⁴⁰ with updated parameters that fit current estimates such as for the Sun Galactocentric radius or the circular velocity curve⁴¹. We used as initial conditions a set of test particles distributed in Galactocentric radius according to a skewed normal distribution with a skewness of 10, location parameter of 4 kpc and scale parameter of 6 kpc. The azimuthal angle was fixed between 0° and 50° to mimic a localized perturbation in the disk. For simplicity, all particles were put at the mid-plane with null vertical velocities. The radial and azimuthal velocities were initialized, respectively, following Gaussian distributions centred at 0 with a dispersion of 40 km s$^{-1}$ and centred at the circular velocity at the particular radius with a dispersion of 30 km s$^{-1}$. The particle orbits were computed for 1 Gyr. To these particles we added particles in the disk that had been integrated for a much longer timescale and are therefore well mixed, to

simulate the stellar populations that have not been perturbed. Of all the particles in the simulation, we used only the ones located in a range of 10° in azimuthal angle at the end of the integration, similarly to the data.

**Model for the horizontal resonances.** The model of Fig. 3c is from a test-particle simulation of orbits integrated in a Galactic potential model including a bar⁴³. The axisymmetric part of the potential was from an existing model⁴⁰. The Galactic bar potential was built using Ferrers ellipsoids⁴⁴ oriented with semi-major axes at 20° from the Sun–Galactic centre line, and with pattern speed set to 50 km s$^{-1}$ kpc$^{-1}$, which corresponds to a period of about 120 Myr. The simulation consisted of 68 million test particles with an initial radial velocity dispersion of 30 km s$^{-1}$ at the solar radius. Their orbits were first integrated in the axisymmetric potential model for 10 Gyr until they were approximately fully phase-mixed. Next, the bar potential was grown in $T_{\text{grow}} = 4$ rotations of the bar. More details on the bar potential, initial conditions and integration procedure are specified elsewhere⁴³. Here we used the final conditions after $T_{\text{grow}}$ (about 500 Myr) and eight additional bar rotations (about 1,000 Myr). From all of the particles in the simulation, we used only the ones located in a range of 10° in azimuthal angle centred on the Sun, similarly to the data.

**Code availability.** We have made use of standard data analysis tools in Python. The codes used to generate the toy models and simulations and to compute the orbital frequencies are available from the corresponding author on reasonable request. The code used to compute the orbits for the potential from ref.⁴² is available at https://github.com/PaulMcMillan-Astro/GalPot. The code used to compute Bayesian distances from parallaxes is available in the TOPCAT platform³⁵.

**Data availability.** The datasets used and analysed for this study are derived from data available in the public Gaia archive (https://gea.esac.esa.int/archive). The Bayesian distances for the Gaia sources with radial velocity³⁷ are available at http://www.astro.lu.se/~paul/GaiaDR2_RV_star_distance.csv.gz. The rest of the relevant datasets and toy models are available from the corresponding author on reasonable request.

29. Luri, X. et al. Gaia data release 2: using Gaia parallaxes. *Astron. Astrophys.* **616**, A9 (2018).

30. Chen, B. et al. Stellar population studies with the SDSS. I. The vertical distribution of stars in the Milky Way. *Astrophys. J.* **553**, 184–197 (2001).

31. Reid, M. J. et al. Trigonometric parallaxes of high mass star forming regions: the structure and kinematics of the Milky Way. *Astrophys. J.* **783**, 130 (2014).

32. Schönrich, R. Galactic rotation and solar motion from stellar kinematics. *Mon. Not. R. Astron. Soc.* **427**, 274–287 (2012).

33. Reid, M. J. & Brunthaler, A. The proper motion of Sagittarius A*. II. The mass of Sagittarius A*. *Astrophys. J.* **616**, 872–884 (2004).

34. Lindegren, L. et al. Gaia data release 2: the astrometric solution. *Astron. Astrophys.* **616**, A2 (2018).

35. Taylor, M. B. TOPCAT & STIL: Starlink table/VOTable processing software. *ASP Conf. Ser.* **347**, 29–33 (2005).

36. Astraatmadja, T. L. & Bailer-Jones, C. A. L. Estimating distances from parallaxes. II. Performance of Bayesian distance estimators on a Gaia-like catalogue. *Astrophys. J.* **832**, 137 (2016).

37. McMillan, P. J. Simple distance estimates for Gaia DR2 stars with radial velocities. *Res. Notes AAS* **2**, 51 (2018).

38. Binney, J. & Tremaine, S. Galactic Dynamics 2nd edn (Princeton Univ. Press, Princeton, 2008).

39. Miyamoto, M. & Nagai, R. Three-dimensional models for the distribution of mass in galaxies. *Publ. Astron. Soc. Jpn* **27**, 533–543 (1975).

40. Allen, C. & Santillan, A. An improved model of the galactic mass distribution for orbit computations. *Rev. Mex. Astron. Astrofis.* **22**, 255–263 (1991).

41. Irrgang, A., Wilcox, B., Tucker, E. & Schiefelbein, L. Milky Way mass models for orbit calculations. *Astron. Astrophys.* **549**, A137 (2013).

42. McMillan, P. J. The mass distribution and gravitational potential of the Milky Way. *Mon. Not. R. Astron. Soc.* **465**, 76–94 (2017).

43. Romero-Gómez, M., Figueras, F., Antoja, T., Abedi, H. & Aguilar, L. The analysis of realistic stellar Gaia mock catalogues – I. Red clump stars as tracers of the central bar. *Mon. Not. R. Astron. Soc.* **447**, 218–233 (2015).

44. Ferrers, N. On the potentials of ellipsoids, ellipsoidal shells, elliptic laminae and elliptic rings of variable densities. *QJ Pure Appl. Math* **14**, 1–22 (1877).

45. Eggen, O. J. Stellar groups. II. The *ζ Herculis*, *ε Indi* and 61 Cygne groups of high-velocity stars. *Mon. Not. R. Astron. Soc.* **118**, 154–160 (1958).

46. Blaauw, A. Remarks on Local Structure and Kinematics. *Symp. IAU* **38**, 199–204 (1970).

47. Skuljan, J., Hearnshaw, J. B. & Cottrell, P. L. Velocity distribution of stars in the solar neighbourhood. *Mon. Not. R. Astron. Soc.* **308**, 731–740 (1999).

48. Antoja, T., Figueras, F., Fernández, D. & Torra, J. Origin and evolution of moving groups. I. Characterization in the observational kinematic-age-metallicity space. *Astron. Astrophys.* **490**, 135–150 (2008).
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Extended Data Fig. 1 | Velocities of the stars at the solar Galactocentric radius.** Two-dimensional histograms of combinations of radial, azimuthal and vertical Galactic cylindrical velocities for the stars in our sample of Gaia data located at 8.24 kpc < R < 8.44 kpc, in bins of 1 km s⁻¹. V_R and V_φ are positive towards the Galactic anticentre and the direction of Galactic rotation, respectively. The darkness is proportional to the number of counts. a, Although the bimodality seen here, separating the Hercules stream from the rest of the distribution, was known⁴⁵,⁴⁶, as well as some other elongated structures in this velocity projection¹¹,⁴⁷,⁴⁸, the numerous and thin arches are a new phenomenon revealed by Gaia data¹². The

semi-circular dotted line marks an arbitrary line of constant kinetic energy in the plane $E_k = (V_R^2 + V_\phi^2)/2$, as predicted for the substructure generated in horizontal phase mixing⁷,⁸. **b**, The data have a box-like appearance, where the extent in $V_Z$ of the arches varies with $V_\phi$ (arrows), probably created by the correlation between the spiral shape and $V_\phi$ seen in Fig. 1c. **c**, Although some velocity asymmetries where noticed before in the $V_\phi-V_Z$ projection¹¹ and attributed to the Galaxy warp, the sharp shell-like features involving $V_Z$, especially at $V_Z \approx -30$ km s⁻¹ and $V_Z \approx 25$ km s⁻¹, were not previously evident. These shells are different projections of the snail shell pattern of Fig. 1a.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 2** | Location of the stars in the sample. **a**, **b**, Two-dimensional histograms with bins of 0.05 kpc in the X–Y (a) and X–Z (b) projections of our sample of Gaia data. The dotted lines mark the selection of stars in the solar Galactic ring between radii of 8.24 kpc and 8.44 kpc. The Sun is located at (X, Y, Z) = (−8.34, 0, 0.027) kpc and the Galactic centre (GC) is marked with a black dot.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 3** | **Modelled vertical positions and velocities of stars with time.** The plots show the snail shells created in the phase space evolution under an anharmonic potential. **a–c**, Phase-space evolution at different times ($t=0, 10, 100, 200, 1,000$ Myr) for an ensemble of particles at a fixed Galactocentric radius of $R=8.5$ kpc with an initial Gaussian distribution in $Z(t=0)$ with mean of –0.1 kpc and dispersion of 0.04 kpc

and in $V_Z(t=0)$ with mean of –2 km s⁻¹ and dispersion of 1 km s⁻¹.
**d–f**, Same as **a–c**, but for a skewed normal distribution of initial radius with skewness of 10, location parameter of 8.4 kpc and scale parameter of 0.2 kpc. In all cases, the evolution is under an anharmonic oscillator derived from the expansion of a Miyamoto–Nagai disk for small Z. In **a** and **d** the stars are colour-coded by vertical period.
<figure><img src="image_12.png" /><figcaption>Extended Data Fig. 4 | Vertical frequency for orbits in a Galaxy model.</figcaption></figure>

a, b, Frequencies as a function of Galactocentric radius *R* computed in the updated model from ref. ⁴¹, colour coded by the vertical amplitude (a) and by the vertical velocity amplitude (b) of the orbits.
<figure><img src="image_13.png" /><figcaption>**Extended Data Fig. 5** | Position of the spiral turns in the vertical positions and velocities. The Z–Vz plane for stars at Galactocentric radii of 8.24 kpc to 8.44 kpc, coloured as a function of median guiding radius Rg in bins of ΔZ = 0.04 kpc and ΔVz = 2 km s⁻¹, with vertical and horizontal lines showing the approximate locations of the observed snail shell (turn-around and mid-plane points).</figcaption></figure>
Extended Data Table 1 | Time estimates from the turn-around points of the spiral

<table><thead><tr><th>Z<br>(kpc)</th><th>ν<br>(rad/Myr)</th><th>Time<br>(Myr)</th></tr></thead><tbody><tr><td>-0.59 ±5</td><td>0.058±0.002</td><td rowspan="2">461<sup>+183</sup><sub>-105</sub></td></tr><tr><td>-0.23 ±5</td><td>0.072±0.002</td></tr><tr><td>0.40 ±5</td><td>0.065±0.003</td><td rowspan="2">566<sup>+290</sup><sub>-140</sub></td></tr><tr><td>0.75 ±5</td><td>0.054±0.001</td></tr></tbody></table>

The first column indicates the vertical positions of the turn-around points, which are equal to the amplitude of the orbits except for the sign, and the estimated uncertainty ranges. The other columns are the frequencies corresponding to these amplitudes and the starting times of the phase-mixing process corresponding to each pair of consecutive spiral turns.
Extended Data Table 2 | Time estimates from the mid-plane points of the spiral

<table><thead><tr><th>A<sub>Vz</sub><br>(km s<sup>-1</sup>)</th><th>v<br>(rad/Myr)</th><th>Time<br>(Myr)</th></tr></thead><tbody><tr><td>-37 ±3</td><td>0.066±0.002</td><td rowspan="2">505<sup>+253</sup><sub>-132</sub></td></tr><tr><td>-19 ±3</td><td>0.079±0.002</td></tr></tbody></table>

The first column indicates the vertical velocities at the mid-plane passages, which are equal to the velocity amplitudes of the orbits except for the sign, and the estimated uncertainty ranges. The other columns are the frequencies corresponding to these amplitudes and the starting time of the phase-mixing process corresponding to the pair of consecutive spiral turns.