Galactic seismology: the evolving “phase spiral" after the Sagittarius dwarf impact

In April 2018, the ESA Gaia astrometric mission (Perryman et al., 2001; Prusti et al., 2016; Soubiran et al., 2018) presented its second data release (DR2). These observations have revolutionized both stellar and Galactic studies (DR2: Gaia Collaboration et al., 2018). After only two years of observations, the remarkable precision of Gaia’s measured stellar parameters has led to new discoveries and new fields of study. These data are being complemented inter alia by large spectroscopic surveys – RAVE (Steinmetz et al., 2006), APOGEE (Majewski et al., 2017), Gaia-ESO (Gilmore et al., 2012), LAMOST (Deng et al., 2012) and GALAH (De Silva et al., 2015). In all respects, this is the golden age of Galactic archaeology (Freeman & Bland-Hawthorn, 2002).

One of the most important discoveries to emerge from the ESA Gaia astrometric survey is an unexpected phase-space signal in the local stellar disc (Antoja et al., 2018). In Galactic cylindrical coordinates ($R,\phi,z$), individual stars have velocities ($V_{R}$, $V_{\phi}$, $V_{z}$) and oscillation frequencies ($\Omega_{R},\Omega_{\phi},\Omega_{z}$) $=$ ($\kappa,\Omega,\nu$). In a volume element defined by ($\Delta R,R\,\Delta\phi,\Delta z$) = ($\pm 0.1$, $\pm 0.1$, $\pm 1$) kpc³ centred on the Sun, the Gaia team detect a coherent spiral pattern¹¹1This feature has been variously referred to as the “snail” (Antoja et al., 2018), the “snail shell” (Li & Shen, 2020; Li, 2020), the “Gaia spiral” (Darling & Widrow, 2019a), the phase-plane spiral (Binney & Schönrich, 2018) and the phase-space spiral (Michtchenko et al., 2019). Consistent with the traditional use of “phase mixing” rather than “phase-space mixing,” we adopt the term “phase spiral” (e.g. Khanna et al., 2019; Hunt et al., 2019; Mackereth et al., 2019; Wang et al., 2019). The advantage of the more compact language becomes apparent when used as an adjective, e.g. phase-spiral evolution, phase-spiral dynamics, etc. in the phase plane defined by $z$ and $V_{z}$. The phase spiral is most evident when each point of the $z-V_{z}$ phase plane is represented by either $\langle V_{R}\rangle$ or $\langle V_{\phi}\rangle$, averaged over the local volume. (Thereafter, Khanna et al. 2019 showed the phase spiral emerges more clearly when encoded by $\langle L_{z}\rangle$, which measures stellar angular momentum about the Galaxy’s spin axis.) This phenomenon is indicative of a system that is settling from a mildly disturbed state to a stationary configuration through the process of phase mixing (Lynden-Bell, 1967).

Since its recent discovery, up to a dozen independent research papers seek to explain or shed light on the phase spiral phenomenon. Explanations to date invoke three different mechanisms: (i) partial coherence imposed by the recent dissolution of several star clusters (Candlish et al., 2014; Michtchenko et al., 2019); (ii) vertical oscillations excited by the strong buckling of a stellar bar (Khoperskov et al., 2019); and (iii) vertical and in-plane oscillations induced by a massive disc-crossing perturber (Binney & Schönrich, 2018; Darling & Widrow, 2019a). The first two mechanisms are of internal origin; the third external mechanism has received the majority of attention to date, not least because the Sgr dwarf galaxy is observed to be undergoing disruption as it crosses the disc. This is the focus of our attention because we consider it to be the most likely interpretation for the phase-spiral phenomenon (Laporte et al., 2019; Bland-Hawthorn et al., 2019). The effect is now observed over too large a surface area (see below) to be a coherent phenomenon arising from dissolving star clusters. Moreover, the extremely tight alignment of the plane of the stellar bar with the inner disc plane (Bland-Hawthorn & Gerhard, 2016) imposes strong constraints on the idea of a buckling bar.

Three independent observing campaigns confirm the phase-space pattern and study it in more detail. With the aid of the GALAH survey, Bland-Hawthorn et al. (2019, hereafter B19) dissected the phase spiral as a function of orbit actions, stellar ages and stellar abundances. They were able to detect the pattern over a larger radial range ($7<R<9$ kpc) and azimuthal angle ($\delta\phi\lesssim\pm 15^{\circ}$). The LAMOST survey provides some evidence for a changing phase-spiral pattern over the radial range $6<R<12$ kpc, as predicted in action-based, analytic models (e.g. B19, their Fig. 20) and observed in the highest contrast numerical simulation to date (Laporte et al., 2019, hereafter L19). Both Li & Shen (2020, their Fig. 4) and B19 present evidence that the phase spiral is more clearly defined in the cold orbits and less clear in the hotter, $\alpha$-enhanced disc, further underscored by the 2-component disc model in B19.

In these studies, the spiral pattern in phase space is found to be somewhat patchy and uneven. This is due in part to the well-documented kinematic (Widrow et al., 2012; Williams et al., 2013; Carlin et al., 2013) and density asymmetries (Yanny & Gardner, 2013; Slater et al., 2014; Xu et al., 2015) when comparing the upper and lower sides of the Galactic plane. However, when the $z-V_{z}$ plane is encoded as the angular momentum about the $z$-axis, $\langle L_{z}\rangle$ (action $\langle J_{\phi}\rangle$), the phase-spiral pattern is smooth and varies systematically as a function of disc location (Khanna et al., 2019). This suggests that the emergence of the phase spiral is somehow related to the angular momentum properties of the disc, in accord with the BS18 model.

Galactic seismology²²2For a brief history of this emerging field, see Appendix A. The terminology is used interchangeably with galactoseismology (Widrow et al., 2012)., a subset of Galactic archaeology, is useful for exploring past interactions and structural properties of the Galaxy. Specifically, we seek a better understanding of how a differentially-rotating stellar disc responds to a powerful impulse and how this sets up the phase-spiral phenomenon. A key question is what are the correct properties to ascribe to Sgr $-$ its total mass, infall velocity, gas fraction and impact radius $-$ at a specific epoch in the past (e.g. Vasiliev & Belokurov, 2020; Vasiliev et al., 2021). In our view, new observations indicate that Sgr’s original mass, presumably prior to infall, is likely to have been comparable to the LMC at the present epoch. In the core of the dwarf, the high metallicity tail in [Fe/H] (Nidever et al., 2020, their Fig. 14) extends well beyond the SMC metal distribution and is comparable to, or exceeds that, of the LMC. The metallicity spread is indicative of the high-mass limit for dwarf galaxies (Lee et al., 2008; Hayes et al., 2020). Its star formation history (Ruiz-Lara et al., 2020) indicates that Sgr lost its gas on successive disc crossings, as modelled for the first time by Tepper-García & Bland-Hawthorn (2018). Sgr’s present total mass is about $2-4\times 10^{8}$ M_⊙ within a radius of 5 kpc from the dwarf’s centre (Law et al., 2005; Vasiliev & Belokurov, 2020). Thus, Sgr’s original total mass and baryon content (Sec. 6) must have been $1-2$ orders of magnitude higher at some time in the distant past (Jiang & Binney, 2000, their Fig. 4).

For consistency with most earlier work, especially Binney & Schönrich (2018), we adopt an initial perturber mass of $2\times 10^{10}$ M_⊙ to drive the impulse (see also Gómez et al., 2013). Earlier work has already shown that this is the low-mass threshold needed to drive the necessary response (L19, B19); models using a lower perturber mass can account for Sgr’s tidal tails (e.g. Law et al., 2005) but produce a very weak vertical disc response. We do not view the apparent inconsistency with observations as a problem. If the phase spiral phenomenon is long-lived, a disc crossing prior to the current event may in fact be responsible. But this raises new questions about the stripping rate experienced by Sgr as a function of time, issues we address in Sec. 6.

In what follows, we start by reviewing the Binney & Schönrich (2018, hereafter BS18) model (Sec. 2), which provides the impetus for our new work, before discussing the overall strategy (Sec. 3). The main goal is to design and carry out a numerical simulation that bridges between the simplicity of an analytic model (BS18) and the complexity of a realistic numerical simulation (e.g. L19). In Sec. 4, we discuss the adopted parameters of the Galactic model and the intruder. The Galactic model is set up carefully with specific initial conditions that guarantee long-term stability in an N-body simulation. In Sec. 5, with reference to web links to our simulations, we describe our analysis and the new results. In the final discussion (Sec. 6), we show that the new simulations are a major improvement over earlier work (L19, B19), although many questions remain.

BS18 develop a purely analytic (toy) model to investigate the response of the Galactic stellar disc to the perturbation induced by an external massive object. They employ action-based modelling defined by distribution functions $f(\mathbf{J})$ that are analytic functions of the action integrals ${\mathbf{J}}=(J_{R},J_{\phi},J_{z})$. With recourse to the Action-based GAlaxy Model Architecture (agama) package (Vasiliev, 2019), they produce a highly self-consistent, multi-component model of the Galaxy. Their vertical density profile through the solar neighbourhood is an excellent match to the observed two-component disc profile (Gilmore & Reid, 1983). This is important if the model is to match the vertical actions ($J_{z}$) and frequencies ($\Omega_{z}\equiv\nu$) given the anharmonic nature of the vertical gravitational potential. They choose phase-space coordinates for a million stars in the solar neighbourhood and backwards integrate them to a time when the perturber first appears. The intruder is moving with a downward velocity of $\vec{v}_{P}\approx 300$ km s^-1, has a mass of $M_{P}=2\times 10^{10}$ M_⊙ and impacts the disc at $R_{P}=18$ kpc along the Sun’s radius vector. They are then able to compute the acceleration imposed by the intruder on the particles relative to the unperturbed acceleration.

BS18 observe that the largest off-planar perturbation is expected at $\phi-\phi_{o}\sim\pm\pi/2$ because these stars are either moving towards ($\phi-\phi_{o}<0$) or away ($\phi-\phi_{o}>0$) from the perturber as it transits the disc. (Here, $\phi_{o}=0$ because the intruder transits along the Sun’s radius vector.) Thus, stars approaching $\phi=\phi_{o}$ feel a net upward pull, and stars receding from $\phi=\phi_{o}$ feel a net downward pull. This action sets up a strong low-order bending mode across the disc. As the perturber transits the disc, stars near the crossing point experience a strong in-plane tug towards the perturber’s trajectory, which binds the in-plane and vertical actions.

With the aid of Figs. 1-3, we illustrate a key insight from the BS18 model. The coordinate system of the $z-V_{z}$ plane is shown in Fig.  1 superposed on the surfaces of section for different stellar orbits confined to tori about the Galactic Centre. In Fig. 2, sloshing in the $z-V_{z}$ plane induced by the impact leads to a stellar overdensity that is offset with respect to the original stratification at a mean angle of $\theta_{z}\approx\theta_{z,o}$ (Fig.  2). Whereas before stars had constant values of $J_{z}$ that were independent of the conjugate angle $\theta_{z}$, the entire distribution now depends on $\theta_{z}$ and the assigned $J_{z}$ associated with the new ellipse. Since for a disc potential, $\Omega_{z}$ declines with $\sqrt{J_{z}}$ (e.g. Candlish et al., 2014, BS18), a phase-spiral structure starts to develop. Crucially, BS18 note that the phase spiral is not apparent in stellar overdensity, although L19 detect a weak signal, in agreement with the original study (Antoja et al., 2018). Therefore, a different mechanism is needed for the pattern to emerge so clearly in $\langle V_{\phi}\rangle$ and $\langle V_{R}\rangle$, which we now discuss.

In Fig. 3, we illustrate how elliptic orbits from the inner and outer disc are able to enter the local volume. The vertical oscillation frequency $\Omega_{z}$ declines smoothly with increasing radius because of the disc’s exponential decline in the local surface density. Stars from the inner disc reach us at their aphelia, and these are characterised by lower $\langle V_{\phi}\rangle$ and higher $\langle\Omega_{z}\rangle$; the converse is true for stars coming in from the outer disc.

In Fig. 2, the dislodged stars now have new positions in the $z-V_{z}$ plane. At every point in that plane, stars have a spread of orbit properties where their clockwise movement in the plane is determined by their current $\Omega_{z}$ value. The dislodged stars trace ellipses in the $z-V_{z}$ plane with a conjugate angle that grows with time $t$ as $\theta_{z}=\Omega_{z}t+\theta_{z,o}$. But $\Omega_{z}$ is tied to $\langle V_{\phi}\rangle$ (and $\langle V_{R}\rangle$) such that stars with lower $\langle V_{\phi}\rangle$ move ahead of stars with higher $\langle V_{\phi}\rangle$. Thus, the phase spiral is defined by a change in $\langle V_{\phi}\rangle$ (and $\langle V_{R}\rangle$) as we move clockwise around each ellipse in the $z-V_{z}$ plane. The change in projected $\langle V_{\phi}\rangle$ (and $\langle V_{R}\rangle$) becomes decoupled from the projected density along the ellipses in Fig. 2. This is why the phase spiral (with its intrinsically higher dynamical range) emerges so clearly when encoded kinematically, and is visible even when the projected density of stars is comparatively low (see Sec. 5.3).

The motivation for our new work can be stated as follows. We seek to understand the phase-spiral phenomenon with the simplest possible model that is true to the BS18 toy model while incorporating the advantages of the N-body approach. First, our work demands a more rigorous treatment for setting up the initial conditions. We require the coldest possible synthetic disc that comes into equilibrium at the start of the simulation and does not develop significant substructure or change its properties over the course of many disc rotations. Darling & Widrow (2019a) demonstrate that the emergence of the phase spiral is a competition between phase mixing and the disc’s self-gravity. Consistent with our approach, a lower self-gravity leads to a stronger signal (higher contrast with respect to the background population) although a higher self-gravity appears to sustain the signature for longer.

Secondly, we require a better treatment of the impulse approximation by an N-body solver. In particular, the heavy point mass must be ‘live’ on its first and only approach, and during the transit, and thereafter it must not continue to disturb the disc. In our earlier work (B19), the disc response can be more erratic in the post-transit phase due to the sustained force applied by the intruder. While there is an in-built asymmetry to our approach, this is true of all Sgr-based orbit families because of dynamical friction and the perturber’s mass loss along its orbit (Law et al., 2005; Nichols et al., 2014). The underlying mechanisms involved in the disc response to an extended, disrupting perturber are obfuscated by many overlapping processes (e.g. L19, B19).

In contrast, the BS18 toy model uses a theoretical description of the Galaxy that is analysed in terms of distribution functions defined in action space, as discussed in Sec. 2. The latter model, while computationally efficient, lacks the realism of the actual system, most notably the complex recoil of both Sgr and the Milky Way to their mutual attraction and, importantly, the self-gravity and ‘elasticity’ of the dominant stellar disc. How the vertical and in-plane motions emerge and work together during the interaction are unclear (D’Onghia et al., 2016; Darling & Widrow, 2019a; Poggio et al., 2020). But by using a model that bridges the divide between analytic and N-body models, we seek to shed more light on this process.

If we are to compare the BS18 model meaningfully with an N-body simulation, the basic parameters of their set-up must be matched closely (Sec. 4.1). In our model, a point mass is used in place of an extended perturber to maintain the comparison and is placed on a hyperbolic orbit to minimize ongoing interactions with the disc. Aguilar & White (1985) note that in order to recreate the impulse approximation within an N-body simulation, the interval during which the mutual gravitational forces are significant is short compared to the crossing time within each system (see also Binney & Tremaine, 2008, §8.2). The “impulsiveness” is further ensured (see Sec. B) by a perturber mass that is constant until the point of disc transit, and then declines exponentially in time such that the mass is negligible at the second disc crossing (Sec. 4.1). In the next section, we explore a new suite of models to simulate a clean impulse imparted by a disc-crossing perturber.

In order to achieve the fidelity of the phase spiral in the BS18 models, these authors recommend that any N-body simulation use $\sim 10^{8}$ disc particles, an order of magnitude larger than any matched, impulsive disc simulation to date. BS18 observe a synthetic phase spiral over an area in the $z-V_{z}$ plane that is comparable to Antoja et al. (2018), i.e. $\pm 60$ km s^-1 and $\pm$1 kpc, a fourfold improvement over L19 and B19. After following BS18’s recommendation, we also achieve phase spiral patterns over the smaller region.

The Gaia discovery of the phase spiral was entirely unanticipated (Antoja et al., 2018). As described in Sec. 4.1, the most natural interpretation is incomplete phase mixing after a major disturbance of the disc. But the remarkable signature raises important questions: What is the precise mechanism that explains it? Is the signature fragile or robust? How widespread is it and can we learn from it?

On the balance of evidence presented, Galactic seismology, a subset of Galactic archaeology, has a rich future. But our work raises more questions than answers. In the closing comments of the last section, we make the case for follow-up studies. Within the broad remit of Galactic seismology, an important goal is to age-date the origin of the phase spiral. Using two distinct mechanisms, B19 and Khoperskov et al. (2019) show that once the corrugation is excited, it can last up to about 1.5 Gyr. But repeated disc crossings from a massive perturber may wipe out the phase coherence imposed by the previous crossing (L19, B19). How are we to understand this?

The orbit period today for the low-mass Sgr remnant is likely to be quite short ($\sim$700 Myr, Ibata & Lewis, 1998). In Sec. 1, we make the case for a high initial mass for Sgr $-$ of order the LMC mass ($\sim 10^{11}$ M_⊙) $-$ for which, due to its orbit within the extended Galactic halo, dynamical friction is likely to be effective. The decay time and mass loss rate require knowledge of Sgr’s initial orbital configuration (Jiang & Binney, 2000). But different orbit models lead to a catastrophic mass loss over a time frame of $7-9$ Gyr such that the observed low mass today is entirely plausible. In external galaxies, stellar streams have been modelled with mass loss rates of order 0.7 dex per orbit loop (Amorisco et al., 2015). Thus, given the mass of Sgr today is too low to excite the phase spiral (Sec. 1), we suggest it was triggered by an earlier crossing in the past 1$-$2 Gyr and that Sgr has been losing mass at a high rate ($\sim 0.5-1$ dex per orbit loop) in the last few passages.

Over timelines of order 1 Gyr, the phase-spiral signal appears to be quite fragile, coming into view in one time step (e.g. Fig. 12) and then either evolving in situ, migrating in azimuth (co-moving RSR frame), or fading altogether in the next time step. In our extended simulation run, we find that the last vestiges of the phase spiral are gone by $T=1.5$ Gyr. Although no panel in Fig. 12 shows a close match in both $\langle V_{\phi}\rangle$ and $\langle V_{R}\rangle$ simultaneously, we find features that are ubiquitous and reminiscent of the discovery paper (see Fig. 15). Note here that the simulated panels show weak evidence of a second phase spiral emerging out of sync with the dominant one. This can arise from having two phase-spiral arms moving through the same volume, which becomes increasingly likely as the density wave becomes highly wrapped at late times.

There is a strong case for repeating the simulations with more particles (N $\gtrsim$ 10⁹) allowing for smaller sample volumes. This would also help improve the simulation’s declining spatial resolution at larger Galactic radius beyond $R\approx 10$ kpc (see Sec. 4.2).

The simulated phase spiral encoded by $\langle V_{R}\rangle$ is acceptable, but the phase spiral encoded by $\langle V_{\phi}\rangle$ shows up to one less wrap than observed. Khanna et al. (2019) found that the Gaia phase spiral is particularly prominent and well behaved in $\langle L_{z}\rangle$, but we were unable to see any improvement in our simulated signal when encoded in the same way. Given the central role of orbital angular momentum in generating the phase spiral (Sec. 2), it is not clear why this should be. While our simulations achieve the recommended particle count of $N\sim 10^{8}$ (Binney & Schönrich, 2018), our test simulations with smaller intrinsic spatial resolution (Sec. 4) did not alter the outcome significantly. But more particles would allow us to better match the Gaia RVS volume and will be necessary for interpreting future Gaia data releases.

There is one important consequence of our “rollercoaster” model that has been glossed over. The interaction of density waves and bending modes implies that the younger stellar populations are those that are more likely to be caught up in the phase-spiral action. It is probable that older stars also contribute to the signal, but the feature is likely to be more enhanced for the younger stars. While there is existing evidence that this is true (B19, Li & Shen, 2020), the dominant population of the phase spiral needs to be addressed more carefully (e.g. better stellar ages) in future studies (Sharma et al., 2016; Miglio et al., 2020).

This brings us full circle to the original prediction of the phase spiral (with two arms rather than the singular arm observed by Gaia) arising from dissolving star clusters (Candlish et al., 2014) that predated the Gaia era. Spiral arms in the Milky Way can be identified from star-forming regions (e.g. Lépine et al., 2001), as they can be in external galaxies, but this is no guarantee that young star clusters ($\lesssim 100$ Myr) are entirely restricted to the spiral arms. An inventory of almost 2000 star clusters confirmed by Gaia (Cantat-Gaudin et al., 2020, Fig. 8) shows that star clusters of all ages inhabit our neighbourhood within a few kiloparsecs, with the youngest star clusters confined to broad bands rather than ridge lines. Most clusters that form appear to dissolve within 100 Myr (Krumholz et al., 2019; Adamo et al., 2020) such that cluster dissolution can conceivably contribute to the innermost wraps of the phase spiral (cf. Li & Shen, 2020). Cluster dissolution is unlikely account for the signature in toto, not least because the models to date predict a 2-arm phase spiral (Candlish et al., 2014), and that has not been observed so far.

Intriguingly, Michtchenko et al. (2019) identify parts of the local phase spiral as arising from several moving groups, including Pleaides, Hyades, Sirius (Ursa Major) and Coma Berenices. These are four of the closest clusters, all of which appear to be unbound. What is particularly striking is that their ages range from 100 to 600 Myr (Famaey et al., 2008) with the youth carrying the stronger signal. More work is needed to assess the detailed structure of the local corrugation (e.g. Friske & Schönrich, 2019; Beane et al., 2019) and its relationship to the local spiral arms (e.g. Miyachi et al., 2019; Griv et al., 2020). How these relate to the local young clusters and star-forming regions also becomes important (Quillen et al., 2020; Cantat-Gaudin et al., 2020).

It is likely that our cold-disc simulations, while very carefully constructed, are too simplistic. Grion Filho et al. (2020) have called for a holistic approach to modelling and understanding galaxy interactions. In this respect, there is a great deal to be done. Vasiliev et al. (2021) have indicated that the LMC’s interaction with the Galaxy and Sgr are important to address, particularly in relation to the development of Sgr’s tidal loops. In the next phase, simulations must introduce a clumpy dark matter halo, giant molecular clouds and smooth gas, a central bar and internally-driven spiral arms, with the overarching goal of producing a more realistic framework for future study. This work will inform the nature of future studies, leading to better targetted surveys rather than all-sky surveys to test different non-equilibrium dynamical models.

To date, we have not considered the contribution of the cool gas. Another prediction from B19 (their Sec. 8.2) is that the corrugation may exist in the Galactic H${\scriptstyle\rm I}$ or the molecular hydrogen distribution, and may even explain Gould’s belt and other local phenomena. Intriguingly, such a gas wave may already be evident (Alves et al., 2020). Connecting this information with star clusters at different stages of their evolution may also provide important insights (e.g. Quillen et al., 2020). If the spiral arms and bending modes are triggered by Sgr, we predict that the corrugations of the wrapping-up wave pattern do not align with the spiral arm pattern (e.g. Fig. 11).

In the spirit of Iye (1985) and Widrow et al. (2012), there is a case for extending galactoseismology to the study of nearby galaxies (cf. Gómez et al., 2020). Arguments have been made for undulations in the Galactic rotation curve (Martinez-Medina et al., 2019) being correlated with the spectacular $R-V_{\phi}$ ridges seen by Gaia (Antoja et al., 2018; Khanna et al., 2019), or possibly even with the spiral arms (Williams et al., 2013). With careful disc modelling and subtraction, it may be possible to identify the signatures of bending modes in substantial numbers of nearby disc galaxies (cf. Matthews & Uson, 2008; Gómez et al., 2020). In such instances, it will be important to establish whether the bending mode wraps up with the disc rotation, as we show here. Morphologically, the wrapping of the spiral density wave is more easily determined and thus it may even be possible to establish whether the interaction of these two wave modes is a regular phenomenon.

The main goal of this work is to provide a better understanding of the response of the Galactic stellar disc to a strong impulse in an N-body setting. Our work shadows Binney & Schönrich (2018) who use an analytic model of a disc-crossing satellite ($2\times 10^{10}$ M_⊙) to explain the phase spiral and its coupled behaviour in $\langle V_{R}\rangle$ and $\langle V_{\phi}\rangle$. Here are the main findings:

1. 1.

   The initialisation code agama (Vasiliev, 2018, 2019) is essential to setting up long-term stability in isolated, multi-component, fully self-consistent realizations (defined at the outset) of the Galaxy, as we demonstrate here. We consider this to be a major step in the field of Galactic seismology.

2. 2.

   As is well known since the advent of N-body codes, a disc-crossing mass drives spiral arms that wrap up as the disc rotates differentially. More recent work has revealed the vertical, wave-like response of the cold disc in addition to the spiral arm perturbation. Here we find that the impulse triggers a superposition of two distinct bisymmetric ($m=2$) modes $-$ a density wave and a bending wave $-$ that wrap up at different rates. Stars in the faster density wave wrap up with time $T$ according to $\phi_{D}(R,T)=(\Omega_{D}(R)+\Omega_{\rm o})\>T$ where $\phi_{D}$ describes the spiral pattern and $\Omega_{D}=\Omega(R)-\kappa(R)/2$. While the pattern speed $\Omega_{\rm o}$ is small, it is non-zero indicating that this is a dynamic density wave. The slower bending wave wraps up according to $\Omega_{B}\approx\Omega_{D}/2$ producing a corrugated wave, a phenomenon now well established in the Galaxy.

3. 3.

   The bunching effect of the density wave triggers the “phase spiral" in the $z-V_{z}$ plane as it rolls up and down on the bending wave (“rollercoaster” model). In other words, the wrapping-up spiral arms undulate with the slower wrap of the bending mode, producing a rollercoaster effect along each spiral arm. This non-equilibrium system is capable of driving the phase-spiral phenomenon.

4. 4.

   In agreement with earlier work (e.g. Antoja et al., 2018), the phase spiral is most evident when each point of the $z-V_{z}$ phase plane is represented by either $\langle V_{R}\rangle$ or $\langle V_{\phi}\rangle$, averaged over the volume, indicating that the disturbance has locked the vertical oscillations with the in-plane epicyclic motions. We find that the phase-spiral action takes about $\Delta T\approx 400$ Myr to become a disc-wide phenomenon, but once it sets in, it can survive till about 1.5 Gyr after the disc transit. (The spiral arms and the bending mode are not strictly bisymmetric with respect to the Galactic Centre: the time delay $\Delta T$ reflects the time it takes the solar neighbourhood and the impact site in the outer disc to realign radially.)

5. 5.

   In earlier work (B19, L19), it was shown that a massive perturber “resets” the disc by wiping out any existing phase-spiral dynamics from an earlier disc transit. Thus, given the mass of Sgr today is too low to excite the phase spiral (Sec. 1), we suggest it was triggered by an earlier crossing in the past 1$-$2 Gyr and that Sgr has been losing mass at a high rate ($\sim 0.5-1$ dex per orbit loop) in the last few passages.

JBH is supported by an ARC Australian Laureate Fellowship (FL140100278) and the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO-3D) through project number CE170100013. JBH & TTG acknowledge the resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. We thank E. Vasiliev for his valuable assistance and continued conversation with the use and development of the agama code. We thank James Binney for earlier conversations that forced us to look at agama carefully. We are indebted to Romain Teyssier for assisting us with the idea of a time-dependent impulsive mass, both in terms of our additions to the code, and for allowing sink particles for the first time to operate with pure N-body simulations. We thank Shourya Khanna and Sanjib Sharma for providing the Gaia DR2 data necessary to recreate the discovery images from Antoja et al. (2018). We are particularly grateful to Ken Freeman for many years of collegiality and for insightful and helpful comments on the manuscript. This work has made used of Matplotlib, a Python-based plotting package (Hunter, 2007), Pynbody (Pontzen et al., 2013) and Gnuplot, originally written by Thomas Williams and Colin Kelley.

The ramses and agama codes were developed by independent investigators and are readily available online. Our work concentrated on designing suitable setup files to run with these facility codes; these are freely available upon request. All simulated data sets specific to the figures are also available upon request.