Snails across Scales: Local and Global Phase-mixing Structures as Probes of the Past and Future Milky Way

Signatures of vertical disequilibrium have been observed across the Milky Way’s (MW’s) disk. These signatures manifest locally as unmixed phase spirals in z–v z space (“snails-in-phase”), and globally as nonzero mean z and v z , wrapping around the disk into physical spirals in the x–y plane (“snails-in-space”). We explore the connection between these local and global spirals through the example of a satellite perturbing a test-particle MW-like disk. We anticipate our results to broadly apply to any vertical perturbation. Using a z–v z asymmetry metric, we demonstrate that in test-particle simulations: (a) multiple local phase-spiral morphologies appear when stars are binned by azimuthal action J ϕ , excited by a single event (in our case, a satellite disk crossing); (b) these distinct phase spirals are traced back to distinct disk locations; and (c) they are excited at distinct times. Thus, local phase spirals offer a global view of the MW’s perturbation history from multiple perspectives. Using a toy model for a Sagittarius (Sgr)–like satellite crossing the disk, we show that the full interaction takes place on timescales comparable to orbital periods of disk stars within R ≲ 10 kpc. Hence such perturbations have widespread influence, which peaks in distinct regions of the disk at different times. This leads us to examine the ongoing MW–Sgr interaction. While Sgr has not yet crossed the disk (currently, z Sgr ≈ −6 kpc, v z,Sgr ≈ 210 km s−1), we demonstrate that the peak of the impact has already passed. Sgr’s pull over the past 150 Myr creates a global v z signature with amplitude ∝ M Sgr, which might be detectable in future spectroscopic surveys.


INTRODUCTION
Our understanding of vertical structure in the Milky Way has drastically evolved over the past decade.Discussions of the disk have traditionally been centered around equilibrium axisymmetric models -planar and fully phase-mixed vertically and with simple periodic perturbations in azimuth.The emerging field of galactoseismology is shifting that focus towards the mild, but significant departures from equilibrium that are increasingly evident in observations.Vertical asymmetries were first pointed out in three distinct data sets by Widrow et al. (2012); Carlin et al. (2013), and Williams et al. (2013).Widrow et al. (2012) found North-South asymmetry in SDSS DR-8 (Aihara et al. 2011) and SEGUE (Yanny et al. 2009).Carlin et al. (2013) used PPMXL (Roeser et al. 2010) and LAMOST (Cui et al. 2012;Zhao et al. 2012) to conclude that stars above and below the midplane exhibit opposite radial motion.Williams et al. (2013) found similar asymmetry around the solar neighborhood in RAVE (Steinmetz et al. 2006).
These local asymmetries in disk motions were shown to be matched by vertical asymmetries of the disk in space, which were traced to several kpc beyond the Sun by Xu et al. (2015).Coincidentally, Price-Whelan et al. (2015) were finding evidence that structures tens of degrees from the plane and at Galactocentric radii of 15-30 kpc, well beyond the traditional limits of the disk (see Newberg et al. 2002, for discovery papers) nevertheless had velocity trends and stellar population properties consistent with disk membership (see also subsequent work that confirms this interpretation Sheffield et al. 2018;Li et al. 2017;Bergemann et al. 2018).These discoveries suggested the local corrugations of the disk were likely part of a global pattern of bending and breathing modes, as first pointed out by Widrow et al. (2014); Gómez et al. (2013).Simulations of LMC and Sgr-like satellites interacting with a Milky-Wayscale galaxy could reproduce the scales of these perturbations, both locally and globally, supporting the plausibility of this interpretation of local and global-scale asymmetries being associated (Laporte et al. 2018a,b).
The reach and high-dimensionality of the Gaia data sets (Brown et al. 2016;Gaia Collaboration et al. 2018a, 2020) allowed clear confirmation of what these earlier studies were hinting at -the existence of global-scale, vertical ripples coursing through our Galactic disk (Gaia Collaboration et al. 2018b).For the first time, Gaia DR-2 enabled the local vertical asymmetries in position and velocity to be dramatically visualized as a clear z-v z phase-spiral1 (Antoja et al. 2018).The richness of the data have inspired analysis and comparison to simulations on both global (e.g.see projections and visualizations in Schönrich & Dehnen 2018;Kawata et al. 2018;Salomon et al. 2020;Poggio et al. 2018a,b;Laporte et al. 2019;Poggio et al. 2020;Eilers et al. 2020;Friske & Schönrich 2019) and local (Antoja et al. 2018;Binney & Schönrich 2018;Darling & Widrow 2019;Laporte et al. 2019;Bland-Hawthorn et al. 2019;Li 2020) scales.In particular, simulations of a Sgr-like satellite impacting a MW-like disk that were developed to fit the pre-Gaia data were shown to contain analogous manifestation of both the local and global signatures of vertical disequilibrium with similar spatial and velocity scales (Laporte et al. 2019).Other works have reinforced these local-global connections with further views of the data, tailored simulations and analytic models (Xu et al. 2020;Bland-Hawthorn & Tepper-Garcia 2020;Widrow et al. 2020;Bennett & Bovy 2021).
While the impact of a satellite provides a natural explanation for the origin of oscillations perpendicular to the Galactic Plane, it is not the only one.Any vertical perturbation, such as bar buckling (Khoperskov et al. 2019), can plausibly do the same (although stellar ages of the phase-spiral suggest this is not the origin in the Milky Way, see Laporte et al. 2019) .Moreover, once bending motions are excited in the disk, self-gravity can launch further disturbances (Darling & Widrow 2019;Bland-Hawthorn & Tepper-Garcia 2020).
Luckily, the combination of local and global responses should provide multiple constraints on the origin of specific features.The coupled evolution of these spirals is largely driven by the simple process known as phase-mixing (see Darling & Widrow 2019, for some limitations of this interpretation).Phase mixing can occur when stars which are initially at random orbital phases are systematically offset by a perturbation to then be at the same orbital phase.If these stars have a range of orbital properties (e.g.frequencies) they will subsequently spread along the orbit and mix in orbital phase.Figure 1 shows the different orbital frequencies as a function of guiding radius, R g (a proxy for distance from the Galactic center; see §3.1 for definition).In the disk x-y plane phase-mixing is driven by the R g -dependence of azimuthal frequencies, Ω φ , while in the z-v z -plane, phase-mixing is driven by the range in vertical frequencies, Ω z , at any fixed R g (see Figure 1b).Phase-mixing can lead to a well-defined spiral for some time, but ultimately this spiral will wind up to the extent that the population once again appears to have randomly distributed phases.
Conceptually, rewinding and interpreting the z-v z and R-φ spirals simultaneously should give us insight into the timing, strength and location of the perturbation that caused them, adding clarity to our understanding of the nature of that perturbation as well as the properties of orbits across the Galactic disk.This paper explores the feasibility of this ultimate goal.Our study uses both test particle and N-body simulations (described in §2) to explore how multiple signatures can be traced back to offer multiple views of a single event (described in §3).The results are used in a first application of the ongoing interaction between Sgr and the MW (described in §4) and future prospects are discussed in §5.

SIMULATIONS & DATA
In this section we describe the simulations and data that we analyze in subsequent sections.Although this paper is based on Gaia eDR-3 data, we interpret those data in terms of simple simulations.After presenting relevant observations in Gaia data, we primarily use test particle simulations which allow us to fully control the orbit of the satellite and the galactic potential.These are described in §2.1.We also show that our results are robust in the presence of self gravity by including a comparison to a self-consistent simulation (in §4), and the self-consistent model is described in §2.2.
All the simulated MW-like disks that we use in this work generically exhibit orbital properties illustrated in Figure 1. Figure 1a shows how φ-rotation frequency Ω φ , vertical epicyclic frequency Ω z , and radial epicyclic frequency Ω R vary with guiding radius, R g (a proxy for distance from the Galactic center; see §3.1 for definition), in an unperturbed simulated disk with stars on near-circular orbits.

Test Particle Simulations
The test particle models are constructed and evolved as described in Hunt et al. (2019), using the galactic dynamics library galpy (Bovy 2015).The initial condition for the disc of massless particles is sampled from galpy's 3D quasiisothermaldf distribution function (adapted from Binney 2010).The distribution function has an initial scale radius R s = R 0 /3, local radial velocity dispersion σ v R = 0.15v c (R 0 ), and local vertical velocity dispersion σ vz = 0.075v c (R 0 ), where R 0 = 8 kpc and v c (R 0 ) = 220 km s −1 .We evolve the disk in galpy's MWPotential2014 for 7 Gyr allowing it to reach equilibrium.
We present three test particle models with varying parameters for the satellite galaxy as summarised in Table 1.In each model, the satellite galaxy is created with galpy's PlummerPotential as a Plummer sphere of mass M sat and scale parameter of 0.8.We calculate the orbit of the satellite by evolving its 'present day' coordinates backwards in MWPotential2014 while using the ChandrasekharDynamicalFrictionForceroutine to take into account dynamical friction.The satellite potential then follows this orbit forward in time with MovingObjectPotential.The three models differ in: (i) values of M sat , (ii) initial phase-space coordinates for the satellite, (iii) time since present day for which the disk is evolved, and (iv) number of test particles in the disk.We do not include a bar or spiral arms in the model potentials, such that the satellite galaxy is the only perturbing influence on the disk.
Model A consists of a 200 million particle disk interacting with a satellite of mass M sat = 2 × 10 10 M initialized with Sagittarius-like phase-space coordinates from Simbad 2 (Wenger et al. 2000;Karachentsev et al. 2004;Ibata et al. 1997;McConnachie 2012) using galpy's orbit.fromname routine.We use this model in §3.2.Although the satellite is the sole perturber to the disk, our conclusions about vertical disequilibrium signatures in §3.4 will apply broadly to generic perturbations (e.g.bar buckling, spiral arms, etc).The mass of the satellite is chosen to be heavier than the remnant mass of Sagittarius (e.g.Vasiliev & Belokurov 2020) in order to generate a strong response, and it is held constant throughout to simplify our model.
Models B and C are used in §4.2, and consist of a 1 billion test particle disk.The orbit of the satellite in both models is initialized with Sagittarius' present day phase-space coordinates which we set to be (x, y, z) = (17.5,2.5, −6.5) kpc, and (v x , v y , v z ) = 237.9,−24.3, 209.0) km s −1 , following Vasiliev & Belokurov (2020).Model B has a satellite of mass M sat = 3 × 10 9 M and the disk has evolved under the influence of the satellite for 150 Myr (ending at the present day), such that the disk experiences only the final passage of Sagittarius.Model C has M sat = 10 10 M and the disk has been evolved for 4 Gyr (ending at the present day).

Self-Consistent Simulations
We also present a self-consistent model for a qualitative comparison to the test particle Models B & C in §4.2.The initial conditions are generated using galactics (Kuijken & Dubinski 1995), using the parameters of Model MWb from Table 2 of Widrow & Dubinski (2005) which result in a disk which remains stable against bar formation over a period of several Gyr (see Widrow & Dubinski 2005, for a thorough analysis of the isolated MW-like disk galaxy).The model contains ∼ 1.12 × 10 9 self gravitating particles, of which ∼ 2.2 × 10 8 are in the disk, ∼ 2.2 × 10 7 are in the bulge and ∼ 8.8 × 10 8 are in the dark halo.
The initial conditions for the satellite are the same as those for Sagittarius in the L2 model of Laporte et al. (2018b), which is composed of two Hernquist spheres (Hernquist 1990).The first represents the dark matter with M 200 = 6 × 10 10 M , c 200 = 28, M h = 8 × 10 10 M and a h = 8 kpc, and the second the stellar component with M * = 6.4 × 10 8 M and a h = 0.85 kpc.
The combined model is evolved using the Bonsai N -body tree code (Bédorf et al. 2012) for 8.3 Gyr with a smoothing length of 50 pc and an opening angle θ 0 = 0.4 radians.The "present day" snapshot in Figure 11c is chosen based on when the satellite is closest to the current coordinates of Sagittarius (Vasiliev & Belokurov 2020) with respect to the Sun, which happens to be at t = 6.88 Gyr.The model will be released alongside a more detailed analysis in Hunt et al. (submitted).

Data
We select stars from Gaia eDR-3 for which 6-D phase space information (parallax, line-of-sight velocity, sky positions, and proper motions) is available.Following Antoja et al. (2018), we require that parallax ω be positive, and that parallax error σ ω be less than 20% (σ ω /ω < 0.2).We use parallax as a proxy for distance (d = 1/ω), and we limit our sample to stars within 7 ≤ R(kpc)≤ 9.This selection contains ∼ 4.6 million stars.Our aim is to explore the origin and evolution of z-v z spirals and understand how these local features relate to the macroscopic vertical ripples in the disk x-y plane.There are three factors to consider in the response to a disk perturbation: (i) phase-mixing around the disk in x-y following the perturbation creating R-φ spirals ( §3.2.1), (ii) phase-mixing in z-v z following the perturbation to form phase-spirals ( §3.2.2), and (iii) the self-consistent disk response causing additional effects.We make the deliberate choice to focus on the combination of the first two phenomenaphase-mixing across dimensions -and defer the addition of the third to future work (see Darling & Widrow 2019, for some cautionary notes on the limitations of our work).In the test particle simulations, the only perturbation which can cause the onset of spirals is the satellite galaxy crossing the disk.This allows us to isolate the time and spatial scales of phase-mixing without the confusion of multiple sources of perturbations.Hence, our conclusions will apply to phase-mixing following any generic vertical perturbation to the disk, but we are missing the effect of self-gravity which can complicate the picture.
We dissect a local sample of Gaia data and an analogous one in the test particle simulations and uncover multiple phase-spirals within each selection ( §3.1).We explore the global context of these phase-spirals in two ways-first by following their evolution backwards through time in the test particle simulations ( §3.2) to understand when and where they were excited, and subsequently by building a toy model to illustrate the spatial and time-scales of the interaction that excited them ( §3.3).We put together our findings in a combined picture of phase-and physical spirals in §3.4.

Dissecting One Phase-Spiral into Multiple
Figure 2a shows the z-v z -plane for Gaia eDR-3 stars in the solar neighborhood, colored by the log number density.A phase-space spiral is visible in this local sample.Figure 2b is the same as Figure 2a, except the color bar represents the fractional overdensity relative to the mean number density (ρ) at each pixel, ∆ ≡ (ρ − ρ)/ρ (following Laporte et al. 2019).
We can dissect this sample further by exploiting the fact that stars which end up within the solar neighborhood today did not always travel together within the same enclosed volume.Hunt et al. (2020) demonstrated how grouping stars around the disk by azimuthal action J φ (equivalent to the z component of the angular momentum L z in an axisymmetric potential) rather than radius R more clearly separates them into sets that have shared histories.(For global disk samples, further grouping by the angle, θ φ conjugate to J φ , rather than physical angle φ can add further clarity to this separation.)The orbits of stars with similar J φ can be characterized by epicyclic oscillations around the same guiding radius, R g , and are hence associated in space.Moreover, since they have the similar R g , they also have similar azimuthal periods, and hence remain associated over time.The guiding radius of a star is calculated by solving R g = J φ /v c (R), where v c (R) is the circular velocity as a function of R. Throughout the paper, we assume a perfectly flat rotation curve with v c (R) ≡ v 0 = 220 kms −1 , and estimate R g ≈ J φ /v 0 .This approximation of a  constant v 0 = 220 kms −1 has been made for simplicity, and therefore the R g values used in our work will not be exact.However, none of our results will be affected by this assumption.Figure 2c is a histogram of R (orange) and R g (navy blue, filled) of the selected Gaia eDR-3 stars and clearly illustrates that although the local sample is limited by physical distance from the Sun, R g allows us to probe a much wider radial range across the disk (∼0-15 kpc, in this case).The "sample" snapshot, 184Myr after the impact when all 5 Rggroups merge into a local volume at φ ≈ −170 • .This figure makes the point that two physical spirals can be identified extending across the disk-one that winds up as we go forward in time, highlighted by the vz color in panels 4(c-e), and another that winds up as we go backward in time, traced by the various simulated Rg groups progressing from panel 4d back through 4a.The time intervals between the panels were chosen to be non-uniform (although symmetric about t − timp = 0) so that the we can visualize the disk response long before/after the impact (t − timp = ±184 Myr), as well as when the satellite is close to the disk (t − timp = ±40 Myr).
In the top panel of Figure 3a, we repeat the same z-v z visualization of Gaia data as Figure 2b, and split this sample into five R g ranges in the bottom panels, between 4-6 kpc, 6-7 kpc, 7-9 kpc, 9-10 kpc, and 10-12 kpc.
The vertical velocities and positions in each panel are scaled by their respective dispersions to adjust the aspect ratio of the phase spiral.We see from the bottom panels that distinct morphologies of the spiral exist at different R g within the local sample (see also Li 2020).
To perform an analogous split in test particle simulations, we select a 30 • azimuthal wedge in the disk between 6 <R(kpc)< 8 at a time chosen simply by virtue of the fact that it exhibits phase spirals with similar wrapping ( 2 wraps), and significant variation across R g , like in Figure 3a.This time happens to be ∼180 Myr after the second passage of the satellite galaxy through the disk, and we will refer to this as the "sample" time.Figure 3b shows rescaled z-v z for simulated stars within the specified local volume in the top panel, and the bottom panels show rescaled z-v z data for the R-limited simulation sample split into five R g bins, between (i) 3-5 kpc, (ii) 5-6 kpc, (iii) 6-8 kpc, (iv) 8-9 kpc, and (v) 9-11 kpc.We use these five R g bins throughout the remainder of §3.As with the Gaia data, each R g bin in the test particle simulation also reveals a different phase-spiral morphology, which otherwise gets obscured in the 6 <R(kpc)< 8 categorization.

Tracing the Evolution of the R-φ and z-v z Spirals by Rewinding the Test Particle Simulations
The fact that the morphology of the z-v z phase spiral within the local volume depends on guiding radius raises the prospect of using this variation to study Galactic history.In order to explore the utility of the varied morphologies, we track particles in the five R g bins introduced in Figure 3b backwards over time.We analyze their projections first in the x-y plane (forming R-φ spirals), and then in z-v z (forming phase-spirals).

Physical Spirals
illustrates the evolution of two types of R-φ spirals, one that winds up as we go forward in time, and another that winds up as we go back in time.The figure shows snapshots starting at 184 Myr prior to the disk passage (panel 4a), through the time of the satellite impact, t imp (panel 4c), to the sample time 184 Myr after the disk passage (panel 4e).The first R-φ spiral is traced by the color across the face of the disk, which represents the mean v z of all the particles.The influence of the satellite from under the disk pulling the particles downwards (blue color, v z < 0) and subsequently upward after crossing the midplane (red color, v z > 0) can be seen in the three middle panels (4b-d).The global response eventually winds up (i.e.phase-mixes) into a clear spiral across the x-y plane by the "sample" time.Note that this simple description of the satellite's impact followed by phase-mixing misses the additional effect present in reality and in the test particle simulation, that the disk particles are also oscillating vertically.We will return to this in the next section.
The colored contours projected onto each panel trace the evolution of the second R-φ spiral.Each contour encloses 50% of the stellar population in one of the five R g groups shown in Figure 3b.In panel 4a, the five groups (starting from the lowest R g group in cyan, increasing through green, purple, magenta, and red), trace a tightly-wound physical spiral.The variation in azimuthal frequencies (see Ω φ (R g ) in Figure 1a) causes this spiral to unwind in each of the subsequent panels until all the groups coincide at the "sample" time in the rightmost panel.It is striking to see in panel 4c that at the time of disk-crossing (t = t imp ), the different groups are spread widely across the disk in azimuth and radius.
We conclude that, because the now-local stellar population was spread across the disk in the past, the z-v z spiral in each local R g group contains distinct information about any past perturbation.We further point out that any local volume in the disk will contain multiple phase-spirals, and the physical spiral in v z across the x-y-plane is a signature of these multiple viewpoints averaged together.

Phase-Spirals: Local Phase-mixing in z-vz
Having tracked the x-y location of the five R g groups over time, we now explore the evolution of their morphologies in the z-v z plane.As noted in §1, while the variation in Ω z with J φ (or R g ) causes a variety of spiral morphologies to be apparent in the same local volume, it is the spread in Ω z at a certain J φ that leads to the phase spiral itself.Figure 1b illustrates both of these points with a 2D histogram of Ω z as a function of J φ for the simulated disk stars.The median Ω z changes with J φ , and at any given J φ there is a range (∆Ω z ) of values present.For instance, at ∼ 10 kpc, the 1σ dispersion about the median σ Ωz ≈ 2 epicyclic orbits per Gyr (marked by the dashed orange curves in Figure 1b).Hence, over a few hundred Myrs we would expect particles with R g ∼10 kpc that are displaced vertically by some perturbation to wind up, with the faster oscillating stars at the head of a z-v z spiral, and the ones with lower Ω z forming the tail.These spirals will fade when particles have had time to fully phase-mix.They will form more slowly in coordinates where stars have smaller spreads in frequencies.
In order to trace the onset and scale of responses in z-v z in the simulated R g groups, we adopt a simple asymmetry parameter (based on that used in Widrow et al. 2012;Bennett & Bovy 2021).The asymmetry (A X ) in a phase-space property X of a group of stars is given by where N (X ≥ 0) (N (X ≤ 0)) is the number of stars which have property X ≥ 0 (X ≤ 0).For our purpose of detecting asymmetry specifically in the z-v z plane, we introduce a combined kinematic asymmetry parameter . Each row ((a)-(e)) in Figure 5 corresponds to a specific time relative to t imp (same time instances as the panels in Figure 4).In each row, the leftmost panel has the time marked, and shows A z,vz as a function of R g ; five panels to the right ((i)-(v)) show the z-v z plane for each of the five simulated R g groups at the respective time.Colored contours enclose 50% of the stellar sample in each R g group.Note that these contours do not necessarily enclose the same stars as the contours in Figure 4. Well before the disk crossing at t−t imp = 184 Myr (row 5a), simulated stars in all R g groups are unperturbed and almost symmetrically distributed in the z-v z plane (i.e., the 50% contours are circles centered on (v z =0, z=0) and A z,vz (R g )≈ 0).By t = t imp (Figure 5c), the 50% contours are distorted and A z,vz has become large (especially for group (v) (red) as it is the one closest to the point of impact; see Figure 4c).We see in rows 5(c-e) that the phase-spirals develop at different rates for the various R g groups.A z,vz oscillates over the course of the disk crossing.We emphasize here that A z,vz is not a measure of how developed (or wound) the phase spirals are, rather, it is a metric for how asymmetrical the distribution of stars is about (v z =0,z=0) .
In Figure 6, we show A vz , A z , and A z,vz as a function of t−t imp for each of the five simulated R g groups.The vertical oscillations of the different groups are captured in Figure 6a, with the period of oscillation being longer for larger R g (as seen in Figure 1c).Figure 6b attempts to capture the growth of the response over time by plotting the amplitude of z-v z distortions overall.As expected, A z,vz ≈ 0 for all R g at early times.We assume A z,vz ∼ O(0.1) indicates the onset of the satellite's influence, which becomes apparent as early as ∼ 100 Myr before t imp .The vertical lines with ' ' symbols mark the time of the maximum kinematic asymmetry for a particular R g group, which corresponds to the time just before the spiral has wound up enough that the asymmetry starts to decline because of phase-mixing.This can occur up to ∼ 50 Myr after t imp , demonstrating the fact that responses can vary significantly with R g .4b), Az,v z begins to grow.The 50% contours are slightly distorted due to the disk's response to the approaching satellite.(c) t = timp: Az,v z grows significantly for group (v) (red) as it is the one nearest to the location where the satellite crosses the disk (see Figure 4c).(d) t−timp= 40 Myr (corresponding to Figure 4d).(e) t−timp= 184 Myr, the "sample" snapshot (corresponding to Figure 4e).In rows (c)-(e), it is apparent that spirals develop at different rates for each Rg group.The asymmetry parameter Az,v z oscillates over time, and is an important indicator of the amplitude of the response, which peaks even before any z-vz spirals appear.Az,v z does not reflect how developed (or wound-up) a z-vz spiral is, but rather is a metric for asymmetric distribution of stars about (vz = 0, z = 0).Note that z and vz for all Rg groups have been rescaled by the dispersions at t−timp= −184 Myr.These are given by σz = {0.41,0.44, 0.46, 0.51, 0.56} kpc and σv z = {44.5,37.4, 34.3, 29.8, 23.4} kms −1 for the five Rg groups respectively.
A striking feature of Figure 6 is the large asymmetry amplitude of the outermost R g group ((v) 9 < R g (kpc) < 11, red curves) compared to the other groups.There are two contributing factors which explain this effect: (1) this R g group is the one closest to the disk crossing region (red, v z > 0) at t =t imp (see Figure 4c), and (2) the vertical oscillation period T z ≈ 100 Myr at R g ≈ 10 kpc (see Figure 1; Ω z /2π(R g = 10 kpc) = 10 Gyr −1 or = 10 epicycles per  1).(a) Av z (t−timp) (solid lines) and Az(t−timp) (dashed lines).Asymmetry oscillates with longer time period for groups with larger Rg.(b) Az,v z (t−timp) for the simulated Rg groups.The colored vertical lines topped with a ' ' mark the time of maximum Az,v z for the Rg groups, and the maxima occur over a span of 50 Myr.These panels reiterate that stars in different regions of the disk experience the same perturbation (a satellite disk-crossing, in the case of our simulations) differently because the amplitude and oscillation frequency of asymmetry varies with Rg.The most important takeaway from this figure is that the effects of the disk-crossing begin to appear ∼ 100 Myr before timp, and last well after timp as well.Even though spirals don't appear until much later, the disk is differentially warped by the satellite's pull for 300 Myr over the entirety of a single passage.
Gyr), is approximately a third of the timescale of the disk crossing (∼ 300 Myr-the time over which the satellite causes significant asymmetry, estimated from Figure 6), thus leading to an enhanced, resonant response.
Overall, we conclude that not only does each R g group experience the interactions from a different location in the disk, but the interactions for each group also occur at different times and with different durations.Moreover, the interaction is far from impulsive, but rather comparable to the orbital times.

Toy Model of the Influence of a Satellite During a Disk Crossing
In this section, we use a toy model of the influence of a satellite crossing the disk in order to place the results of prior sections in context -how the experience of the same satellite perturbation can vary across R g groups.For the purpose of isolating the scale of the influence in different regions of the disk, we assess the impact on toy model stars moving in the midplane (z fixed to 0) on perfectly circular orbits (v fixed to φ) throughout the encounter and ignore their vertical oscillations or any vertical displacement due to the satellite's pull.
Figure 7a sketches the toy model: a galactic disk (shown face-on) comprises stars on clockwise circular orbits in the midplane, and a satellite on a vertical trajectory passing through it with v z > 0 (out of the page).We can get some intuition for our results by first considering the three toy model stars shown in Figure 7a, one with φ (t imp ) = φ sat at the time of disk crossing (yellow star), a second with φ (t imp ) < φ sat (blue star), and the third with φ (t imp ) > φ sat (red star).Figure 7b demonstrates that the yellow star experiences exactly equal and opposite force from the satellite before and after t imp , and thus experiences net ∆v z = 0.The blue star is closer to the satellite when it is being pulled down (at t < t imp ), but farther when being pulled up (at t > t imp ) and thus has net ∆v z < 0. Finally, the red star is farther from the satellite when it is being pulled down (at t < t imp ), and closer when being pulled up (at t > t imp ) and thus has net ∆v z > 0. Figure 7c simply reiterates the asymmetric response by showing the cumulative a z (t − t imp ) which is 0 at t − t imp = t i (t i t imp ) for the three sample toy model stars, and by t − t imp = t f (t f t imp ), cumulative a z is positive for the red star (φ sat − φ (t imp ) > 0), 0 for the yellow star (φ sat − φ (t imp ) = 0), and negative for the blue star (φ sat − φ (t imp ) < 0).
We quantify the scales of the satellite influence across the disk by integrating the z-acceleration it exerts on some sample particles over the course of the encounter.Each test particle has a position vector in cylindrical coordinates, r = (R , φ (t), z = 0), such that its galactocentric radius is constant and it remains in the midplane.The azimuth of a particle is given by φ = φ(t imp ) + v t.The satellite (pink circle) of mass M sat is on a vertical trajectory with constant galactocentric radius, azimuth, and upward vertical velocity, described by r sat = (R sat , φ sat , z sat = v z,sat t)).The vertical acceleration a z (t) of a particle due to the gravitational pull of the satellite is given by, (2) Note that we are intentionally only considering the acceleration due to the satellite (and not including the restoring force of the disk) because we want to isolate the response generated by the satellite.Thus in this paper, a z refers to the vertical acceleration defined in eq(2), and we neglect the evolution of orbits by forcing the toy model stars to remain on circular orbits.We calculate a z (t) for all particles in a toy disk where we set M sat = 2 × 10 10 M , z sat (t) = v z,sat t, v z,sat = 339 kms −1 , φ (t) = φ (t imp ) + v t, and v (R ) is obtained from the rotation curve of MWPotential2014 in galpy.The satellite's mass, vertical velocity, and disk crossing coordinates are chosen to closely match the satellite in Model A of the test particle simulation at t =t imp .Over a certain time period t i ≤ t ≤ t f , we can find the total change in a toy model star's vertical velocity, (3) Figure 8a shows a toy disk colored by ∆v z integrated over −0.5 < t − t imp (Gyr) < 0.5.The x-y coordinates of stars in this panel are frozen at t imp to show net ∆v z at t = t f as a function of position at t = t imp with respect to the satellite at disk-crossing.Toy model stars with φ sat − φ (t imp ) < 0 experience net ∆v z < 0 over the entire course of the disk crossing and lie in the blue region, whereas stars with φ sat − φ (t imp ) > 0 have the opposite response and lie in the red region.Only at φ sat − φ (t imp ) = (0, π) is ∆v z = 0. Three sample radii of R= (4, 8, 16) kpc with four sample toy model stars at each radius are chosen (shown with ' ' symbols on dashed circles in Figure 8a) to examine a z (t) for different radii and azimuths.Each of Figures 8b-d shows a z (t − t imp ) for each of the sample radii, and the colored curves corresponds to a ' ' in Figure 8a of the same color at that radius.The horizontal bars mark T φ and T z for the corresponding R. Note that epicyclic oscillations are not included in the toy model; T z is simply taken from the simulated data in Figure 1a to compare typical oscillation timescales with the disk crossing timescale.
The toy model captures the qualitative diversity in the disk response to a satellite perturbation.It demonstrates the R-and φ -dependence of a z (t).It underscores the fact that the satellite's integrated influence on vertical motions: (i) can have opposite signs in different disk regions; (ii) reaches its peak at different times for different azimuths at the same radius; (ii) extends over timescales comparable to both the orbital and vertical oscillation times.Of course, these conclusions are a function of the particular passage we choose to examine -these are the characteristics we expect for satellites with disk crossings at radii that are comparable to the disk size itself.Another caveat of this simple model is that it neglects ∆v R and ∆v φ of stars, which will qualitatively alter the signatures of vertical disturbance we see in the toy disk.

Combined Implications of Our Results
We have explored the origin of the multiple z-v z spiral morphologies apparent when local samples in the Gaia data and test particle simulations are divided into R g bins.In Section 3.2.1 we showed how the locations of these distinct samples, traced back in time, outline an R-φ spiral across the disk (see distribution of colored contours in Figure 4a).Furthermore, the various z-v z spiral morphologies in any local volume are local signatures of a global v z spiral spanning the extent of the disk (see v z structure in Figure 4e).
In Section 3.3 we used a toy model to examine the case of a 2 × 10 10 M satellite crossing the Galactic disk plane and demonstrated how the overall influence is neither local nor impulsive.Stars at different R g which form phase spirals "today", started responding to a perturbation tens of Myrs before t imp and well before coherent phase spiral structures developed.The combination of these results imply that the multiple local phase-spirals, even if associated with the same event, cannot be simply "rewound" to a single time and location to learn about the impact.Rather, each R g group represents distinct viewpoints of the same event which are widely spread out in space and time.
We reiterate that although we present the example of a disk-and-satellite interaction as the cause of these phasespirals, their subsequent phase-mixing in z-v z and R-φ depends only on the disk properties and not the nature of the interaction.Distinct morphologies in local z-v z samples would occur with other types of perturbations as well, and contain multiple viewpoints on the cause, whether a buckling bar, rippling disk or an impacting satellite.Hence local spirals could be powerful diagnostic tools of global disk disturbances.The orbit has been integrated in MWPotential2014 from galpy.The current position of Sgr is marked by the large red dot, the dashed curve is the the past orbit, and the dotted curve is the future orbit.The black dots are 50 Myr apart.The black horizontal line represents the MW disk, with the Sun marked by the orange ' '.The asymmetry of Sgr's past and future orbit indicates that soon after crossing the midplane, its distance from the disk will increase significantly compared to the past ∼ 150 Myr while it was sweeping under the disk.Its proximity to the disk in the recent past leads us to explore whether signatures of vertical disturbance have already developed.

RESULTS II: APPLICATION TO THE ONGOING MW-SGR INTERACTION
The results from previous sections show that a Sagittarius-like satellite passing through the midplane at R≈ 15 kpc can influence the disk at times ∼ ±150Myr around t imp .Currently, Sagittarius is approaching the MW from z ≈ −6 kpc with v z ≈ 200 kms −1 and is expected to hit the outer disk at R ∼ 18 kpc in ∼30Myr.This suggests that, even though Sagittarius has not yet crossed the midplane, there could be signatures of this encounter already developing in the disk.
Figure 9 shows the recent and near future path of Sagittarius's by tracing its orbit within galpy's MWPotential2014 forwards and backwards from its present-day Galactocentric phase-space coordinates from Vasiliev & Belokurov (2020) noted in §2.1.The geometry of Sagittarius's orbit causes it to travel quite close under the disk plane as it approaches its present location rather than simply passing vertically.To explore possible signatures of the current interaction, we now analyze global deviations from zero in v z as well as local phase-space signatures which might be detectable in future surveys.On these short timescales (∼150 Myr), we expect disk self-gravity to be least important and our toy and test particle models to capture much of the response.

Estimates of Scale from the Toy Model
We apply our toy model to explore the nature of Sagittarius's past, present, and future influence on the MW throughout the current disk passage.The toy disk once again comprises particles on strictly circular orbits rotating clockwise with constant circular velocity v (R ) and z = 0, but this time for a satellite travelling along the orbit whose x − −z projection is shown in Figure 9.We set M Sgr to a constant 3 × 10 9 M (the initial Sagittarius mass in Vasiliev & Belokurov 2020).
Figure 10a shows the disk x-y plane at present day (t = t 0 , "today"), colored by ∆v z (see eq(3)) integrated over −0.5 ≤ t − t 0 (Gyr) ≤ 0. The black dashed curve shows the x-y projection of Sagittarius's past orbit, with the black dots spaced at intervals of 50 Myr.The large red dot marks Sagittarius's current x-y position, the red dotted curve is its future orbit and the red cross is where it will cross z = 0.The orange ' ' at (−8, 0) marks the Sun.The plot shows that particles across the disk receive a net velocity kick of up to several kms −1 from this portion of the passage alone.Again, remember that the toy model neglects epicyclic oscillations and thus is not a prediction for the mean v z , but rather for the scale of Sagittarius's influence in different regions.Moreover, the mass of Sagittarius that we have used (M Sgr ≈ 3 × 10 9 M ) is the initial and largest mass in Vasiliev & Belokurov (2020) (see their Figure 9).The Sgr remnant loses significant mass by present day, which would mean that the strength of the v z signal in reality will be much weaker than the v z amplitude in our toy model.
Figures 10b-d each show a z (t) for radii R=4, 8, and 16 kpc respectively and particles at four different azimuths (' ' symbols in 10a).The asymmetry of a z (t) around t 0 emphasizes the fact that most of the vertical perturbation within R 10 kpc due to Sagittarius's ongoing interaction with the Milky Way has already happened over the past 200Myr, as the satellite sweeps close under the disk.More specifically, the ratio between |∆v z | induced over the past 200 Myr and |∆v z | induced over t 0 ± 200 Myr , averaged over the 4 sample toy model stars at each radius is 0.82, 0.79, and 0.69 for R = 4, 8, 16 kpc respectively.

Morphological Predictions from Test Particle Simulations
Our qualitative conclusions from the toy model in §4.1 motivate searching for a signature of the vertical response to the ongoing MW-Sagittarius interaction.We use test particle simulations of a disk galaxy perturbed by a Sgr-like satellite on the orbit prescribed by MWPotential2014 with present day Sgr phase-space coordinates from Vasiliev & Belokurov (2020).We analyze the simulations for two different cases-Model B: a MW-like disk which has only evolved under Sgr's influence over the past 150 Myr with M Sgr = 3 × 10 9 M , and Model C: a MW-like disk that has experienced multiple crossings of Sagittarius over the past 4 Gyr, with M Sgr = 10 10 M .This helps us extricate the vertical response over the past 150 Myr from the remnant effects of past disk crossings.It also allows us to see how the amplitude of the response scales with M Sgr .

Physical Spirals
In Figures 11(a,b), we show a comparison of the x-y plane at present day colored by v z ×10 10 M /M Sgr for Models B and C.There is a distinct pattern in v z which appears to be more or less consistent between the two panels, indicating that it is largely a result of the ongoing MW-Sgr interaction, and not a remnant of past disk crossings.The maximum v z amplitude induced in the disk is significantly lower than the toy model (| max v z ,toy | × 10 10 M /M Sgr ≈ 7/0.3 kms −1 ≈ 23 kms −1 ) because epicyclic oscillations average out the signal.Moreover, the amplitude of v z ∝ M Sgr , such that maximum |v z | × 10 10 M /M Sgr ≈ 3 kms −1 remains almost the same in both panels.2) for purely illustrative purposes to show that a similar pattern in vz emerges even when self-gravity of the disk is accounted for.Note that the x-y limits of this plot are different from panels (a) and (b), which means that the extent of the vz pattern is not the same.However, this can be explained by the various differences between the Sagittarius orbit and disk potential between the test particle and self-consistent simulations.The main takeaway from this figure is that in all three panels, a blue region with vz < 0 emerges around φ ∼ 0 • , and a red region with vz > 0 is present around φ ∼ −90 • .This leads to a robust prediction that these signals in vz are present in the MW disk today and might be detectable in future surveys.Finally, Figure 11c shows a self-consistent simulation snapshot (see §2.2 for details of the model) at present day which has evolved over the past 6.88 Gyr, and the current Sgr remnant mass is M Sgr,rem ≈ 8 × 10 9 M .Although the x-y extent of Figure 11c is different from panels 11(a,b), in all three models, a blue (v z < 0) patch is visible around φ ∼ 0 • on the opposite side of the disk across from the Sun ((x , y ) = (−8, 0)), and a red (v z > 0) region emerges with the highest amplitude of positive v z around φ ∼ −90 • .These two regions are roughly marked by black contours in panel 11a, and appear around the same azimuth in panels 11(b,c) as well (albeit in different x-y positions).The blue patch (labeled '(1)v z < 0' in panel 11a) is created by the downward pull of the Sgr-like satellite, it's present-day x-y location marked by a red dot in each panel.The red region (labeled '(2)v z > 0' in panel 11a) consists of simulated stars which were pulled downwards by Sgr ∼ 100 Myrs ago, and are now traveling upwards as part of their vertical epicyclic motion.
The reason why the coordinates and v z amplitude of the blue and red patches differ between the panels 11(a,b) versus panel 11c is that the disk potential and satellite orbit in the self-consistent simulation are different from those in the test particle simulation, and there is little control over these in the self-consistent model.We emphasize that the existence of the v z < 0 and v z > 0 regions not just in the test particle models but also in a self-consistent MW-Sgr interaction is main point of Figure 11.The fact that a self-consistent disk also exhibits a significant v z < 0 signal near the location of the dwarf is a strong indicator that this signature might exist in our own Galaxy.If detected in future surveys (like the SDSS-V Milky Way Mapper; Kollmeier et al. 2019), the amplitude could be used to infer Sagittarius's remnant mass precisely while the shape could be used to infer trends in disk frequencies (and hence the force field) in that region.There is the caveat that the remnant mass (M Sgr ≈ 8 × 10 9 M ) inducing the v z signal in Figure 11c is much larger than the present day M Sgr quoted in Vasiliev & Belokurov (2020).Thus it is possible that Sgr's tidal effects are negligible at present compared to the self-sustained bending waves in the disk.

Phase-Spirals
We do not find phase spirals in the test particle disk models due to the ongoing Sagittarius-MW interaction, as might be anticipated since there is little time for these features to develop.However, there is significant vertical asymmetry in the simulated disk caused by the interaction, which we present in Figure 12.Panels 12(a,b) respectively show the z-v z plane for particles in the blue region (labeled '(1)v z < 0' in panel 11a) and the red region (labeled '(2)v z > 0' in panel 11a).The black contours in 12(a,b) enclose 50% of the simulated sample in each of the two regions.In panel 12a, the sample is clearly offset toward v z < 0, and a significant fraction of simulated stars have z < 0. Panel 12b shows the black contour offset toward v z > 0 (as expected) and z < 0, indicating that although stars in this patch have started to travel upward with positive v z , they are still below the midplane.

CONCLUSIONS
The main conclusions of this paper are the following: 1.As shown in Li (2020), selection of stars based on their azimuthal action J φ (or equivalently, their guiding radius R g ) makes it possible to resolve multiple z-v z spiral morphologies (Figure 3).We make the case that each of these local phase-spirals probes distinct regions of the disk and began developing at distinct times (Figure 5).These multiple z-v z spirals can originate from the same perturbative event (Figure 4), and therefore each of them offers a different perspective on the same event.
Multiple z-v z spiral morphologies (over a wide range in R g ) exist in the local sample because the effects of a single satellite disk-crossing are long-lasting (∼ 300 Myr) and affect the entire disk.The varied spirals are a reflection both of the fact that different regions of the disk experience the perturbation with different amplitudes (demonstrated with a toy model in Figure 8) and that different regions respond with different characteristic frequencies.Coincidence of the disk-crossing timescale and orbital time period can further amplify the distortions (Figures 1,6).
Since a single perturbative event can have such an extended impact on disk dynamics, we investigate the ongoing MW-Sagittarius interaction.Our investigation leads to several insights about the current disk passage.
2. Even though Sgr is not expected to cross the midplane for another ∼ 30 Myr, the bulk of the influence on the inner disk (R 10 kpc) from this imminent passage has already happened (Figure 10).
3. We do not find z-v z spirals in test particle simulations of the ongoing interaction (Figure 12), but significant asymmetry is expected in R-φ, leading to a disk-wide physical spiral (seen in Figure 11).The amplitude of this v z signature scales linearly with M Sgr .
4. The fact that there is a v z < 0 (blue) patch around Sgr's "present-day" x-y position, and a similarly sized v z > 0 (red) patch around φ ≈ −90 • in both the test particle and self-consistent disk (Figure 11), suggests it is likely that these signatures are present in our Galaxy as well.These patches might be detectable in future Milky Way surveys (e.g.SDSS-V Milky Way Mapper; Kollmeier et al. 2019).If the true mass and orbit of Sgr are indeed such that there is a coherent v z offset, the amplitude can be used to infer M Sgr , while the morphology and location will help constrain both the properties of the disk and Sgr's orbit.However, it is possible that the projection shown in Figure 11c is in reality obscured by pre-existing or independent disk dynamics.
We consider the above conclusions to be generic consequences of phase-mixing alone.A full interpretation of observed features will need to take account of the disk self-gravity, and signatures of the more recent interactions will need to be disentangled from prior perturbative events.Nevertheless, our results demonstrate intuitive starting points towards building methods that can tease apart these overlapping responses.

Figure 1 .
Figure 1.(a) Best-fit curves for epicyclic oscillation frequencies expressed as Ω/2π(Gyr −1 ) for simulated stars in an unperturbed disk in the z, R, and φ directions as a function of Rg=J φ /v0.The Ω(Rg) curves in this figure are fit to near-circular orbits, i.e., for which J 2 R + J 2 z /J φ 1, where JR and Jz are orbital actions in the R and z directions respectively.(b) A 2D histogram showing the spread in Ωz as a function of Rg not restricted to near-circular orbits anymore.The solid orange curve in this panel (different from the orange curve in panel (a)) is the median Ωz, and the dashed curves mark the 1σ dispersion about the median.The spread in Ωz at a particular Rg value, Rg , is what causes a phase-spiral to develop in the z-vz plane.Stars at the head of a z-vz spiral have high Ωz(Rg ), whereas stars in the spiral tail have lower Ωz(Rg ).The variation of Ωz with Rg is what causes multiple spiral morphologies to develop once the disk is perturbed.

Figure 2 .
Figure 2. (a) The phase-space spiral as seen in the z-vz plane with data from Gaia eDR-3.We select ∼ 4.6 million stars within 7 <R(kpc)< 9 with parallax error σω < 20%.vz and z have been rescaled by their respective dispersions, (σv z , σz) = (25 kms −1 ,0.37 kpc).The color represents the log number density.(b) Same as panel (a) except the color shows the filtered number density of Gaia eDR-3 stars, ∆ ≡ (ρ − ρ)/ρ.Note that ∆ is used specifically to highlight the spiral morphology, and it makes the spiral appear as a stream-like structure in phase-space.However, panel (a) is a more accurate representation of how the stars are actually distributed.(c) A histogram of R and Rg (≈ J φ /220 kms −1 ) in our Gaia eDR-3 selection demonstrates that the R-limited sample spans a much wider range in Rg.
(a) Top panel: z-vz data of the local volume of stars taken from Gaia eDR-3 (7 < R(kpc)< 9, same as Figure 2b).Bottom row: Gaia eDR-3 stars in the top panel now split into 5 Rg groups ranging 4 < Rg(kpc)< 12.In all panels, we have rescaled z and vz by the respective dispersions (σv z , σz) = (25 kms −1 , 0.37kpc) to adjust the aspect ratio of the spirals.It becomes clear that a Rg categorization resolves distinct phase spirals (each panel in the bottom row shows a different spiral morphology) which otherwise get averaged out in the R selection (top panel).(b) Top panel: z-vz data of a local sample of stars (6 < R(kpc)< 8) in Model A of our test particle simulation.Bottom row: Simulated stars in the top panel now split into 5 Rg groups ranging 3 kpc< Rg < 11 kpc.In all panels, we have rescaled z and vz by the respective dispersions (σv z , σz) = (27.6 kms −1 , 0.36 kpc) to adjust the aspect ratio of the spirals.Once again, as is the case with the Gaia phasespiral, Rg categorization (bottom row) resolves distinct phase spirals which are averaged out in the test particle simulation R selection (top panel).Colors outlining the Rg ranges specified in each panel (cyan, green, purple, magenta, red) and lower-case roman numerals ((i)-(v)) are identifiers used in figures throughout §3 to refer to the various Rg groups.

Figure 4 .
Figure 4.The 5 simulated Rg groups (shown in Figure 3b) are represented here by contours enclosing 50% of the stars in each Rg bin, overplotted on the disk x-y plane colored by vz.Each panel shows the positions and velocities of simulated stars at a different time around the time of "impact", t = timp (i.e., the time when the satellite crosses the disk midplane).(a) The x-y plane 184 Myr before timp shows how different Rg groups are spread out across the disk.(b) 40 Myr before impact.(c) at the time of impact, t = timp.(d) 40 Myr after the impact.(e)The "sample" snapshot, 184Myr after the impact when all 5 Rggroups merge into a local volume at φ ≈ −170 • .This figure makes the point that two physical spirals can be identified extending across the disk-one that winds up as we go forward in time, highlighted by the vz color in panels 4(c-e), and another that winds up as we go backward in time, traced by the various simulated Rg groups progressing from panel 4d back through 4a.The time intervals between the panels were chosen to be non-uniform (although symmetric about t − timp = 0) so that the we can visualize the disk response long before/after the impact (t − timp = ±184 Myr), as well as when the satellite is close to the disk (t − timp = ±40 Myr).

Figure 5 .
Figure 5.Each row (a)-(e) corresponds to a specific time relative to timp(indicated in the leftmost panels), with the same time intervals as between the panels of Figure 4.The leftmost panel in each row shows the z-vz asymmetry parameter, Az,v z ≡ A 2 z + A 2 vz (explained in the discussion around eq(1)) as a function of Rg.The five plots to the right ((i)-(v)) in each row show the five simulated Rg groups (introduced in Figure 3b) in the z-vzplane.The colored contours in the z-vz panels enclose 50% of the stellar population in each Rg group, and the z-vz asymmetry Az,v z of each group is marked by a dot of the corresponding color and lower case roman numeral in the leftmost panel.(a) Az,v z and the simulated Rg groups at t−timp= −184 Myr.Az,v z is minimal at this time.Most of the stars in all groups are symmetrically distributed in z-vz, as is evident from the 50% contours being circles centered on (vz = 0, z = 0).This row corresponds to Figure 4a.(b) At t − timp = −40 Myr (corresponding to Figure4b), Az,v z begins to grow.The 50% contours are slightly distorted due to the disk's response to the approaching satellite.(c) t = timp: Az,v z grows significantly for group (v) (red) as it is the one nearest to the location where the satellite crosses the disk (see Figure4c).(d) t−timp= 40 Myr (corresponding to Figure4d).(e) t−timp= 184 Myr, the "sample" snapshot (corresponding to Figure4e).In rows (c)-(e), it is apparent that spirals develop at different rates for each Rg group.The asymmetry parameter Az,v z oscillates over time, and is an important indicator of the amplitude of the response, which peaks even before any z-vz spirals appear.Az,v z does not reflect how developed (or wound-up) a z-vz spiral is, but rather is a metric for asymmetric distribution of stars about (vz = 0, z = 0).Note that z and vz for all Rg groups have been rescaled by the dispersions at t−timp= −184 Myr.These are given by σz = {0.41,0.44, 0.46, 0.51, 0.56} kpc and σv z = {44.5,37.4, 34.3, 29.8, 23.4} kms −1 for the five Rg groups respectively.

Figure 6 .
Figure6.Asymmetry parameters as a function of t−timp for the five simulated Rg groups.For details about the different asymmetry parameters shown, see discussion around eq(1).(a) Av z (t−timp) (solid lines) and Az(t−timp) (dashed lines).Asymmetry oscillates with longer time period for groups with larger Rg.(b) Az,v z (t−timp) for the simulated Rg groups.The colored vertical lines topped with a ' ' mark the time of maximum Az,v z for the Rg groups, and the maxima occur over a span of 50 Myr.These panels reiterate that stars in different regions of the disk experience the same perturbation (a satellite disk-crossing, in the case of our simulations) differently because the amplitude and oscillation frequency of asymmetry varies with Rg.The most important takeaway from this figure is that the effects of the disk-crossing begin to appear ∼ 100 Myr before timp, and last well after timp as well.Even though spirals don't appear until much later, the disk is differentially warped by the satellite's pull for 300 Myr over the entirety of a single passage.

Figure 7 .
Figure 7. (a) The x-y projection of a toy disk with a satellite (pink circle at rsat = (Rsat, φsat, vz,satt)) passing through on a strictly vertical path with vz,sat > 0 (out of the page).The disk rotates clock-wise, with all stars remaining strictly in the midplane and at constant radius represented by r = (R , φ (t), z = 0).(b) Vertical acceleration caused solely by the satellite, az(t − timp), where timp is the time the satellite crosses the midplane, for the 3 toy model stars marked in panel (a).(c) Cumulative az(t − timp) for the 3 sample toy model stars.Panels (b,c) together show that over an extended time ti < t < t f such that the initial time ti timp and final time t f timp , the red star (φ (timp) > φsat) gets a net positive ∆vz(t f ), the blue star (φ (timp) < φsat) gets a net negative ∆vz(t f ), and the yellow star (φ (timp) = φsat) experiences a net zero change in vz at t = t f .See discussion in §3.3 for further explanation.(Note: The time axes in panels (b) and (c) are different; we have zoomed into a shorter time interval in panel (b) compared to (c) to show the az(t) curves in better detail.)

Figure 8 .
Figure 8.(a) A toy disk rotating clockwise in the x-y plane colored by vertical velocity ∆vz over −184 < t − timp (Myr)< 184.The positions are frozen at ximp, yimp at t = timp.The white '×' marks the point where the satellite crosses the midplane with vz,sat > 0, coming out of the page.This panel demonstrates the point that toy model stars with φ (timp) > φsat get a net positive kick in vz (red region), the toy model stars with φ (timp) < φsat get a net negative kick (blue region), whereas toy model stars with φ (timp) = φsat or φ (timp) = φsat − π get a net 0 change in vz.For stars in this toy model, ∆vz is calculated by only accounting for az caused by the satellite passing through the disk (see Eq.(2)).That is, we neglect epicyclic oscillations and self gravity.Three sample radii are chosen (black dotted circles at 4, 8, and 16 kpc) with sample stars (marked with ' ' symbols) distributed in azimuth to analyze az(t).(b)-(d) For each sample radius ((b) 4 kpc, (c) 8 kpc, and (d) 16 kpc), we show az(t − timp) for the sample stars selected at that radius.The color of each curve in these panels corresponds to a colored ' ' marked in panel 8a.The red curves have a net positive ∆vz, yellow and green curves have net zero, and blue curves have net negative ∆vz.The black horizontal bars indicate time periods of oscillations at the sample radius (taken from Figure 1a): the longer bar is T φ (period of φ-rotation), the shorter one is Tz (period of vertical epicyclic oscillation).The vertical acceleration of the toy disk stars last over a period comparable to orbital time scales.Note: The stars in this toy model are not oscillating, the horizontal bars are simply shown to indicate what the oscillation time periods are in a MW-like disk at particular R values, and how they compare to the time scale of the disk crossing.

Figure 9 .
Figure9.The x-z projection of Sagittarius's orbit with present day phase space information fromVasiliev & Belokurov (2020).The orbit has been integrated in MWPotential2014 from galpy.The current position of Sgr is marked by the large red dot, the dashed curve is the the past orbit, and the dotted curve is the future orbit.The black dots are 50 Myr apart.The black horizontal line represents the MW disk, with the Sun marked by the orange ' '.The asymmetry of Sgr's past and future orbit indicates that soon after crossing the midplane, its distance from the disk will increase significantly compared to the past ∼ 150 Myr while it was sweeping under the disk.Its proximity to the disk in the recent past leads us to explore whether signatures of vertical disturbance have already developed.

Figure 10 .
Figure 10.This figure's panels and color schemes are the same as Figure 8; MSgr = 3 × 10 9 M .(a) The x-y plane colored by ∆vz integrated over the past 500 Myr up until present day.Since Sgr has been below the disk for the past ∼ 200Myr and is at z ≈ −6 kpc at present, there is a net negative ∆vz across the entire disk.The red dot shows the current x-y position of Sgr, and the black dashed line tracks its past trajectory, with the black dots being 50Myr apart in time.The red dotted line tracks Sgr's future trajectory and the red cross marks the position where Sgr will cross z = 0 (∼30 Myr in the future).The orange ' ' at (−8, 0) kpc indicates the current position of the Sun.(b)-(d) az(t) at different radii and azimuths, with t = 0 at present day.The dotted vertical line marks the time in the future when Sgr will cross z = 0.The extreme asymmetry of az around t = 0 in these three panels makes the point that a major fraction of the ∆vz due to the imminent passage of Sgr has already been induced by present day.Once again, this toy model does not account for epicyclic oscillations and self-gravity, therefore should not be interpreted as quantitatively accurate.

Figure 11
Figure 11.(a) Model B: a test particle disk x-y plane at present day (t = t0) colored by vertical velocity scaled by MSgr vz × 10 10 M /MSgr (this rescaling helps enhance the vz signal from a satellite with MSgr< 10 10 M ).The test particle disk has evolved under the influence of Sagittarius only over the past 150 Myr with MSgr = 3 × 10 9 M .The outlined patches correspond to the panels of Figure 12.(b) the same plot as (a), but for Model C. The test particle disk has evolved over the past 4 Gyr, and MSgr = 10 10 M .The vz pattern across the disk is largely the same between panels (a) and (b), and the vz amplitudes are also similar once scaled by MSgr.The fact that the isolated disk in (a) shows the same vz pattern as (b) implies that the current MW-Sgr interaction is largely what creates this vz signature.(c)"Present day" snapshot from a self-consistent simulation (described in §2.2) for purely illustrative purposes to show that a similar pattern in vz emerges even when self-gravity of the disk is accounted for.Note that the x-y limits of this plot are different from panels (a) and (b), which means that the extent of the vz pattern is not the same.However, this can be explained by the various differences between the Sagittarius orbit and disk potential between the test particle and self-consistent simulations.The main takeaway from this figure is that in all three panels, a blue region with vz < 0 emerges around φ ∼ 0 • , and a red region with vz > 0 is present around φ ∼ −90 • .This leads to a robust prediction that these signals in vz are present in the MW disk today and might be detectable in future surveys.

Figure 12 .
Figure 12.The regions with largest |∆vz| in simulations of the ongoing MW-Sgr interaction do not exhibit phase-spirals.(a) Corresponds to vz < 0 patch outlined in Figure11a.The black contour here encloses 50% of the simulated stellar population in the patch.The sample is clearly offset toward vz < 0, and we also see that a significant fraction of stars in that region have z < 0. (b) Corresponds to the vz > 0 patch in Figure11a.Once again, 50% of the simulated sample is enclosed by the black contour here, which is clearly shifted toward vz > 0 (as expected) and z < 0. This figure additionally provides z information which is not apparent from Figure11a.

Table 1 .
The different models used for test particle simulations (see §2.1 for discussion)