Tidally induced spiral arm wraps encoded in phase space

Context. The impact of Sagittarius and other satellite galaxies such as the Large Magellanic Cloud on our Galaxy and in particular its disc is gradually being disclosed. Simulations tailored to the interaction of the Milky Way (MW) and Sagittarius show rings and spiral arms appearing in the Galaxy disc. However, spiral arms can also be induced by the bar or by disc instabilities. Aims. We aim to study the dynamics of tidally induced spiral arms in the context of the di ﬀ erent encounters with Sagittarius and determine their kinematic signatures in the shape of ridges and waves in angular momentum, similar to those detected with Gaia DR2. Methods. We built toy models of the interaction between a host and a satellite galaxy using orbital integrations after a tidal distant encounter. We derived analytically the shape of the structures seen in phase space as a function of time for simple power-law potential models. We compared these models to a more realistic N -body simulation of the MW Sagittarius-like interaction and also to real data from Gaia DR3. Results. As previously found, an impulsive distant tidal approach of a galactic satellite generates a kick in velocities that leads to a two-armed spiral structure. The arms are made of orbits in between their apocentres and pericentres, thus, they correspond to regions with average negative galactocentric radial velocity. The two-arm pattern rotates at an angular speed of ω − 1 / 2 κ which depends on Galactocentric radius, thus causing winding with time range of 0.8–2.1 Gyr, respectively. This winding produces ridges in the R − V φ projection with alternating signs of V R and oscillations of V R in the L Z − φ space, similar to those observed in the Gaia data. The frequency of these kinematic features increases with time, o ﬀ ering a powerful means to infer the potential and the perturbation’s onset time and azimuthal phase. Fourier analysis allows us to date the impact times of simple models and even to date perturbations induced from subsequent pericentric passages that appear as simultaneous waves. For the MW, the Fourier analysis indicates a superposition of two di ﬀ erent frequencies, conﬁrming previous studies. Assuming that both are due to impulsive and distant pericentre passages, we ﬁnd perturbation times < 0.6Gyr and in the range of 0.8–2.1Gyr. The latter is compatible with a previous pericentre of Sagittarius and would be associated to about four wraps of the spiral arms in the observed radial range. Conclusions. Further work on the self-gravitating response of galactic discs and possible degeneracies with secular processes induced by the bar is necessary. Our study is a ﬁrst step towards shedding more light on the elusive structure and dynamics of the spiral arms of the Galaxy.


Introduction
The hierarchical formation of galaxies in the paradigm of cold dark matter has consequences, not only on how galaxies grow, but also on how they evolve dynamically (e.g. Dressler 1980;Barnes & Hernquist 1992). Undeniable evidence of this are the many galaxies that have been caught interacting with a companion (Arp 1966). The Milky Way (MW) is no exception: the impact of Sagittarius and other galactic satellites such as the Large Magellanic Cloud (LMC) on our Galaxy and, in particular, on its disc is gradually unfolding. For example, host-satellite interactions have been proposed as a mechanism of excitation of the Galactic warp (e.g. Hunter & Toomre 1969;Weinberg 1995;Ibata & Razoumov 1998) and disc heating (Quinn et al. 1993). The vertical asymmetries in the density and in the vertical velocity seen in Widrow et al. (2012) also raised suspicion of disc perturbations from satellites (Gómez et al. 2013). More recent observational evidence for this process came from the discov-ery of the phase spiral  in the second data release (DR2; Gaia Collaboration et al. 2018a) of the Gaia mission (Gaia Collaboration et al. 2016).
Mergers can alter the morphology of galactic discs, exciting the formation of bridge and tail arms in the outermost regions, and of bars, rings, and spiral structure in the inner parts (e.g. Toomre & Toomre 1972;Noguchi 1987;Barnes & Hernquist 1992;Kazantzidis et al. 2008;Pettitt & Wadsley 2018). Simulations tailored to the interaction of the MW and Sagittarius also show that rings and spiral arms would appear in our disc as a response to the interaction (Younger et al. 2008;Purcell et al. 2011). On the other hand, spiral arms can also be excited internally by bars or can be self-excited disc modes.
There are multiple theories about the dynamics, formation, and persistence of spiral structure (e.g. density wave model, swing-amplification theory, transient structures, invariant manifolds; see reviews by Sellwood &Masters 2021 andBaba 2014). For the particular case of the MW, there is evi-dence of spiral arms in young and old populations (Georgelin & Georgelin 1976;Reid et al. 2009;Drimmel & Spergel 2001), but a precise global map of the spiral structure is missing, let alone a solid theory for its origin and dynamics. Persistent unknown aspects are the number of spiral arms in the Galaxy, their exact composition (gaseous versus stellar), their long-or short-lived nature, their excitation mechanism and relation with the bar, and how they rotate with respect to the stellar disc (as a rigid pattern, with pattern speed depending on radius, or co-rotating with the stars; Roca-Fàbrega et al. 2013). Despite all these unknowns, there are fewer doubts about the importance of spiral structure in the dynamical evolution of the stellar disc, its gas, and the newly formed stars (e.g. Debattista et al. 2006;Sellwood 2014).
Density structures such as spiral arms and rings, independently of their origin or dynamics, must be intrinsically related to patterns in the planar velocities. The mean velocities in the X-Y projection around the Sun seen first in the Gaia DR2 data by Gaia Collaboration et al. (2018b) may be related to spiral arms and, interestingly, different types of arms could lead to different velocity patterns (e.g. Grand et al. 2015;Antoja et al. 2016). Some of the observed diagonal ridges in the V φ -R projection Kawata et al. 2018) and EDR3 data (Gaia Collaboration et al. 2021) have been tentatively related to transient spiral structure (Hunt et al. 2018(Hunt et al. , 2019Khanna et al. 2019), to resonances of density wave arms (Hunt et al. 2019; Barros et al. 2020;Michtchenko et al. 2018), or to recent crossings with the spiral arms (Quillen et al. 2018). Other ridges are very likely related to the bar Hunt et al. 2018;Fragkoudi et al. 2019;Monari et al. 2019;Fragkoudi et al. 2020;Laporte et al. 2020). Ridges can also appear due to an interaction with an external satellite (e.g. Laporte et al. 2019;Khanna et al. 2019).
Most likely, different processes are at play in the dynamics of the MW disc. Ramos et al. (2018) found different slopes in the ridges seen in the V φ -R projection, which could indicate different origins. Friske & Schönrich (2019) discovered a complex wave of V R as a function of angular momentum and speculate that it could be a superposition of the effects of the bar and the spiral structure. It is, therefore, of the utmost importance to understand the dynamics of different types of spiral arms and identify observable signatures that help us distinguish them from each other and from other disturbances. In particular, the dynamics of tidally induced arms in the context of the MW remains poorly explored.
Here we study the effects of a perturbing satellite (similar to Sagittarius) on the MW disc, paying special attention to how these structures are seen in phase space in the shape of ridges and other kinematic disturbances similar to those detected in the data. Our aim is to study their dynamics, distinguish them from other types of arms, and ultimately constrain the spiral structure of the MW and the interaction with Sagittarius. We first use simple toy models of impulsive and distant encounters to explore the formation and dynamics of spiral structure (Sect. 2). We link the spiral arms generated in these models to the ridges in the V φ -R projection and to a wave of V R as a function of angular momentum (Sect. 3). Since we see that the frequency of these features becomes larger with time due to the phase-wrapping, we develop a method to recover the time of impact and the parameters of the potential through Fourier analysis (Sect. 4). We then look at the more realistic N-body simulation by Laporte et al. (2018) that we compare against our simpler models (Sect. 5). Finally, we compare models and real MW data from Gaia DR3 (Gaia Collaboration et al. 2022) in Sect. 6. We discuss our results in Sect. 7 and conclude in Sect. 8. Table 1. Models used in this work. The columns indicate: the label of the model, the type of model, the potential used, the type of initial conditions, and the radial velocity dispersion at R 0 in units of V 0 .

Toy models of kinematic waves
In this section we present a series of models of the Galactic spiral arms after an impulsive and distant tidal encounter, in particular, those often referred to as kinematic spiral waves. To introduce basic equations and give a brief historical overview, we start from the idea of the rotated ellipses (Sect. 2.1) and continue with models from orbit integrations for different distribution functions (Sect. 2.2). These are all two-dimensional models (2D) that we compare to the more complex and three-dimensional (3D) N-body model of Sect. 5. In Table 1 we list all these models and specify their details. We use cylindrical coordinates (R, φ, V R , V φ ). In all our models, the plots are presented with rotation in the clockwise direction and φ and V φ are positive in this sense. The origin of φ is set at the negative X axis.
For convenience, in the toy models presented below we use spherical power-law models that in the 2D idealisation are equivalent to an axisymmetric potential with orbits confined to the plane. In these potentials, the circular velocity curve is: where n is the slope of the curve, and V 0 its value at R = R 0 . We explore models with a perfectly flat circular velocity curve (logarithmic potential) and slightly decreasing (n = −0.1) or increasing (n = 0.1) circular velocity curves. The frequencies of nearly circular orbits (with guiding radius R g ) for models of the type of Eq. 1 are: The circular velocity curves as well as the azimuthal ω, radial κ, and ω− 1 2 κ frequencies of these models are shown in Fig. 1. For completeness we add other slopes of the circular velocity curves corresponding to n=0.5, -0.5 (Keplerian), and 1 (homogeneous sphere). As default values we use V 0 = 240 km s −1 and R 0 = 8 kpc, similar to the MW case (e.g. Gravity Collaboration et al. 2021;Reid & Brunthaler 2020;Schönrich et al. 2010).

Lindblad and Kalnajs ellipses
The concept of the kinematic spiral waves was born around the 1950's, in the words of Lindblad (1956), as a "line of thought which may lead to a satisfactory dynamical theory of the spiral phenomenon". The basis of the idea is that, in a reference frame that moves at frequency ω− 1 2 κ, orbits are seen as closed ovals, or equivalently, that the apsidal lines move at this frequency relative to each other if this frequency does not remain the same at all radii (this is the precessing rate of the closed ovals).
To illustrate this, we present our Model 0 in Fig. 2. In this model we built a series of orbits using the equations of perfect epicyclic orbits: For simplicity, in this example we used a flat circular velocity curve (n = 0). In these equations A is the epicyclic amplitude and t is the time. The constant ψ 0 gives the phase in the epicycle but is irrelevant in this particular exercise and it is set to 0. This is similar for φ 0 which is the initial azimuth of the orbits and we took it to be 0. We sampled different R g and times t, choosing A = 0.1, and printed the orbits in the ω − 1 2 κ reference system where they become closed orbits (Fig. 2, top left). These ellipses will appear in this configuration only for a given instant of time as the different ω − 1 2 κ at different radius would make the ovals precess at different rhythms. Lindblad (1956), Lindblad & Lindblad (1958) and Kalnajs (1973) noticed, however, that this frequency can be nearly constant at intermediate and outer radius of disc galaxies (Fig. 1). In this way, the set of ovals remains quite fixed relative to each other creating barred shapes, or, if the ovals appear rotated, persistent trailing spiral structure (left panel of the second row, which has been obtained by rotating the ellipses of the top panel giving them different values of φ 0 , in particular with the inner and outer ellipses differing by 150 deg). We refer to Binney & Tremaine (2008, Sect. 6.2) for more details.
In a more realistic case where ω − 1 2 κ is not exactly constant with radius but has a small slope, a configuration like that of the top panels evolves towards more tightly wound spirality, such as in the second and third rows of Fig. 2 (the latter done by increasing the difference between inner and outer ellipses to 300 deg). This structure slowly curls with time (at a rate of ω − 1 2 κ), which is slower for raising circular velocity curves (right bottom panel of Fig. 1), with the limit of perfect flatness (thus, no curling with time) for the homogeneous sphere (n = 1).
The second and third columns of Fig 2 indicate in colours the radial and azimuthal velocity (the later represented with respect to the circular velocity V c ), respectively, which follow from Eqs. 3. As Kalnajs (1973) already noticed, the spiral arms form by the accumulation of orbits in-between their apocentres and pericentres, thus coinciding with regions where orbits are travelling inwards (negative V R ) and have V φ − V c = 0. Immediately past the arms, the orbits reach their pericentre and have, thus, maximum V φ and V R = 0. The radial and azimuthal velocities show a π/4 rad shift in azimuth. With time, the regions of accumulations of orbits become more concentrated in configuration space, or, in other words, the angular region occupied by the spiral overdensitiies (and consequently the regions of V R < 0) is smaller compared to the inter-arm regions.
Interestingly, a configuration like the one in Fig. 2 (top) can result after a tidal interaction with a companion galaxy. For example, following Struck et al. (2011), the gain in velocities for an impulsive and weak interaction is: where D is a scale parameter (the radius of the disc), and ∆V is a velocity amplitude obtained from the tidal constants and the ratio of the disc size to the distance of closest approach. Spiral patterns of the kinematic-wave type that roll up at a rate ω − 1 2 κ have been reported in multiple N-body simulations (Oh et al. 2008;Struck et al. 2011;Oh et al. 2015;Hu & Sijacki 2018;Pettitt et al. 2016;Pettitt & Wadsley 2018;Semczuk et al. 2017;Hu & Sijacki 2018). This idea has also been recently put in context of the MW by Bland-Hawthorn & Tepper-García (2021) who measured the windup frequency of the spiral structure formed in an N-body simulation following an impulsive impact of a Sagittarius-like perturber. Deviations from the ω − 1 2 κ rate are seen when a bar is also present and can depend on the self gravity of the disc and on time (Pettitt & Wadsley 2018;Oh et al. 2015). In those cases, the frequency usually lies between ω − 1 2 κ and ω − 1 4 κ. We see in Sect. 5 that some of the pericentres of the N-body simulation that we explore also follow these kind of spiral patterns that slowly wind up with time.

Model 1: Full orbital integration
In our modelling series that we dub Model 1, we used galpy 1 (Bovy 2015) to integrate orbits. This is more realistic than our previous Model 0 in which all orbits followed perfect epicycles. We built four different models of this type (Table 1) combining different initial conditions and potentials. We used two different sets of initial conditions from the distribution functions of Dehnen (1999) with radial velocity dispersion at R 0 = 8 kpc of 0.05 and 0.2 times V 0 , which for V 0 = 240 km s −1 give 12 and 48 km s −1 , respectively. The labels of the models (Table 1) indicate the dispersion with c and h standing for cold and hot, respectively. All of them have 200 000 test particles. We sampled the  Kalnajs (1973). The top panel shows a series of orbits in the reference frame of ω − 1 2 κ for which the orbits close. In the first, second and third column, we plot the orbits, the orbits coloured as a function of radial velocity (V R ), and as a function of azimuthal velocity (V φ − V c ), respectively. The two bottom rows shows the same but for increasingly rotated ellipses. distribution only between 0.6 and 1.9 R 0 (4.8 and 15.2 kpc, respectively) in order to focus on the central parts of the disc. We use a scale-length of 1/4R 0 , and an exponential radial-velocitydispersion profile with scale-length R 0 . We chose three analytical potentials from the power law family with n = 0, n = −0.1 and n = 0.1 (marked with solid thicker lines in Fig. 1). In the model labels this slope is indicated with f, r, d indicating flat, rising and decreasing circular velocity curves, respectively. Although we generated distribution functions consistent with the corresponding slope, the initial conditions were integrated for 10 Gyr in order to diminish the effect of the phase-mixing patterns that appear as a consequence of not having a fully consistent distribution function.
After this, we added the velocity kicks of Eq. 4 keeping positions unchanged. We used arbitrary values of ∆V = 10 km s −1 and D = 20 kpc. These directly give the initial conditions for the orbital integration and are represented in the top row of Fig. 3 in the configuration space X-Y. The first column gives the positions of the toy orbits. Second and third columns indicate the average radial velocity V R and azimuthal V φ −V φ (R), respectively, where V φ (R) is the average azimuthal velocity at that radius. The quadrupole signature of Eqs. 4 after the simulated satellite impact at t = 0 is seen, with V R and V φ shifted by a phase of π/4. With this scheme the orbits are distributed at special epicyclic phases depending on their location inside the disc.
The next rows of Fig. 3 show the evolution of the positions and velocities of Model 1hf (with hot initial conditions and the logarithmic potential). The initial pattern rotates at an angular speed of ω − 1 2 κ, which depends on radius R, and thus its locus follows: Fig. 3. Tidal spiral arms of Model 1hf. We show the time evolution between 0 and 1 Gyr (indicated in the panels of the first column) of the disc in the X-Y space, in density (left column), coloured by average radial velocity V R (middle) and average azimuthal velocity V φ −V φ (R) (i.e. with respect to the mean at the given radius, third column). The dashed and dashed-dotted lines show the prediction of the spiral arms location following the winding with ω − 1 2 κ. The sector of the disc marked with dotted lines is a region whose dynamics we examine later in Fig. 6.
Because of the radial dependence of the rotation speed, the pattern winds up (shears) with time. In the top panel of Fig. 3 we show a dashed and a dot-dashed line corresponding to opposite azimuthal regions that have negative radial velocities at the initial times. In the following rows, these lines are wound up as described in Eq. 5. These lines mark the position of stars that are at an epicyclic phase of ψ = π/2 at each time. A material arm (made of the same stars at all times) would instead phase-mix at a rate of ω, that is, with a much steeper dependence on R and thus much faster. Another fundamental difference with material arms is that here the pattern is composed of stars with common epicycle phase at a given time, and thus stars that instantaneously are on top of the crest of the wave change with time.
Continuing with Fig. 3, the epicyclic oscillations lead to regions of high orbital crowding, creating a trailing two-armed spiral structure. In particular, the location of the maximum density (locus of the arms) corresponds to regions with stars whose epicyclic phases passed apocentre and now move towards pericentre. This leads to the correspondence between the locus of the arms and regions with negative radial velocity. This also results in the regions of negative V R getting compressed in configuration space (i.e. their azimuthal width at a fixed radius gets smaller). The inter-arm regions correspond to positive radial velocities. This is the same behaviour as in the formalism presented in Sect. 2.1 (Fig. 2). For this particular model, even after 1 Gyr, Fig. 4. Spiral arms of different models in polar coordinates. We show the φ-R projection at 1 Gyr, in density (left column), coloured by average radial velocity V R (middle) and average V φ − V φ (R) (right). The black lines show the predicted location of the spiral arms.
clear spiral structure is seen in the range of radii from 5-10 kpc (similar to that of the volume sampled by Gaia).
A polar representation often does a better job at showing the spiral arms in density. Figure 4 shows the positions and velocities in the φ-R plane now for different models at 1 Gyr. Again, the azimuthal component of the velocity corresponds broadly to a phase shift of π/4 with respect to V R . The first two rows show different levels of kinematic temperature of the initial conditions from cold to hot (1cf, 1hf). The potential used was the same and, therefore, the winding rate and the exact location of the arms and associated velocity regions are equal. The differences are related to the level of dispersion (noise) in the average positions and velocities.
We note that additional substructure appears both in density and velocities in Fig. 4. For example, a clear pattern of small substructure in density, V R and V φ is seen especially for the cold model (first row). These patterns must be the result of the phase mixing but their exact mathematical descriptions were not derived here where we focus on the structure of our interest (i.e. the spiral arms). We note that similar substructure is observed in similar models such as from Struck et al. (2011).
The three bottom rows of Fig. 4 compare different slopes of the circular velocity curve for the hotter initial conditions: Model 1hf with n = 0 flat (second row), Model 1hd with n = −0.1 decreasing circular velocity curve (third row) and Model 1hr with n = 0.1 rising circular velocity curve (last row). The different models have different slopes of the frequencies with radius, in particular of ω − 1 2 κ: as shown in Fig. 1 this frequency is flatter for the rising circular velocity curve, and becomes progressively steeper for the flat and the decreasing curves. This results in a more tightly wound pattern in the 1hd case (smaller pitch angle).
The degree of winding with time of these kinematic-wave toy models depends, therefore, both on the time (after impact) and on the slope of the circular velocity curve. For a certain time, the number of spiral wraps crossing at a certain azimuth φ i within the range between R min and R max is given by where φ(R) is the one defined by Eq. 5. Figure 5 shows the number of spiral wraps at certain radius ranges as a function of time and of the slope of the circular velocity curve. Since this number depends on φ i , for the figure we show the maximum number of Fig. 5. Number of spiral wraps. We show in colours the expected maximum number of spiral wraps at different radius ranges (4-10 and 10-15 kpc in the first and second columns, respectively) as a function of time t and slope of the circular velocity curve n. We have used 21 for the upper limit of the colour bar but the number of wraps goes beyond this in the bottom part of the left panel.
wraps encountered for all possible φ i . At inner radius, the number of wraps is much larger at a given t and n, as expected from the slope of the curve ω − 1 2 κ. At a certain time and radius, the number of wraps is larger the more negative the slope of the circular velocity curve n is. In the outer parts of this idealised disc, the number of wraps is small even for larger times after the supposed impact.
Below we explore in detail these kinematic waves seen in other projections of phase space. This helps us to link them with the ridges in the R-V φ plane and a V R wave with angular momentum, which are features observed in the MW.

Ridges and angular momentum waves
In Fig. 6 we focus on a particular azimuthal region of the disc to look at other projections of phase space. We start with Model 1hf (hot initial conditions and flat circular velocity curve). We take only the region of φ = [π − 0.2, π + 0.2] rad (marked with dotted lines in Fig. 3), which corresponds to a sector of ≈22.9 deg centred on the positive X axis. In the first column we show a density histogram of the radius of the particles. We see some clear bumps that correspond to locations of higher density, that is to say the spiral arms generated in the model. The vertical dashed black lines show the predicted positions of the arms. These are obtained by first inverting Eq. 5 and using Eq. 2c to obtain: Then we find the positions of the multiple wraps, R m (m = 1, 2, etc) by adding π to the azimuth φ an integer number of times m, that is changing φ to φ + mπ to obtain R m = R(φ + mπ, t). For this particular example we take φ 0 = π/2 (the initial azimuth of the arms in our model) and the azimuth φ = π (the centre of the sector chosen). As shown before, the number of wraps of spiral arms crossing the region becomes larger with time. The spatial separation between wraps at a given time increases with R due to the smaller slope of ω − 1 2 κ with R at large radii. The second column of Fig. 6 shows the density in the R-V φ space, while the third one is the same but coloured by average V R . At the beginning (t = 0, first row), there is a flat distribution of V φ as a function of R because in this model the circular velocity curve is flat. The dispersion in V φ is larger for inner radii, as expected for typical disc galaxies and exactly as generated in the initial conditions. The φ = π ± 0.2 rad region corresponds to initially positive radial velocities (thus dominated by red colours Fig. 6. Ridges and waves for Model 1hf. We show the time evolution between 0 and 3 Gyr (indicated in the first column) of the disc region at φ = π ± 0.2 rad (marked with dotted lines in Fig. 3). The first column shows a normalised density histogram of the radius of the particles R. The second and third columns show the R-V φ projection in density and coloured by mean V R , respectively. The fourth column shows the particles in the L Z -V R projection. The vertical dashed lines in all panels indicate the prediction of the spiral arm locations following a winding rate of ω − 1 2 κ.
in the third column) as can be seen also in the top middle panel of Fig. 3. As time goes by, sawtooth-shaped oscillations appear in the R-V φ projection, with the parts of positive slope roughly corresponding to negative radial velocities (blue colours) and the decreasing parts of the oscillations having positive radial velocities (red colours). As we have seen before, the negative V R regions are linked to the kinematic spiral arms and thus we see that in this projection the spiral arms correspond to the blue bands. We also plot vertical dashed black lines in the radius where the arms are located, which indeed coincide with the centres of the blue bands. The blue and red velocity bands tend to get stepper with time, with the blue ones becoming almost vertical due to the spiral arms becoming more tightly wound (smaller pitch angle). The bands of alternating negative and positive radial velocity also become closer due to the winding of the pattern. Finally, in the fourth column of Fig. 6, we show the radial velocity as a function of the angular momentum. The initial distribution shows simply a straight line with a small positive slope, resulting from the imposed initial perturbation that increases linearly with R -combined with the transformation from R to L Z for a perfectly flat circular velocity curve in an azimuth where particles have not gained or lost angular momentum (panel in the third column of the first row in Fig. 3). Again we see larger velocity dispersion at inner radii. With time, oscillations in V R appear, with increasingly larger frequency. The growing amplitude with L Z comes from the imposed velocities at the beginning (Eq. 4). The black dashed vertical lines correspond to the location of the arms in the space of L Z (corresponding to the minima in V R ). They were obtained similarly as R m . We first transformed Eq. 6 to angular momentum (taking into account that the arms are made of orbits exactly at their guiding radius) and then set φ 0 = π/2, φ = π, and changed φ to φ + mπ, that is computing L Z,m = L Z (φ + mπ, t). For a flat circular velocity curve (n = 0), Eq. 7 becomes: At a given time, the frequency of the V R oscillation depends on L Z , as it also depends on R. This can be seen in Eq. 8 since the separation between valleys of minimum V R is ∆L Z ≡ L Z,m − Article number, page 6 of 19 T. Antoja et al.: Tidally induced spiral arm wraps encoded in phase space L Z,m+1 ∝ 1 m 2 +k 1 m+k 2 , where k 1 and k 2 are constants and m takes larger values for the valleys at smaller L Z (and similarly for n 0 with Eq. 7). However, from Eq. 7 one also realises that the oscillations in V R have constant period if plotted as a function of L −1 Z instead of L Z . We note a similar behaviour for non-flat circular velocity curves as a function of L n−1 n+1 Z . This is explored further in Sect. 4. Figure 7 shows the ridges and waves for the different flavours of Model 1 at 1 Gyr of evolution. In the first two rows we compare the different initial conditions but use the same logarithmic potential model. The effects of the velocity dispersion are observed between panels. Taking the second, third and fourth rows we can compare different slopes of the circular velocity curves (n = 0, −0.1, 0.1). The effect of n is to modify the overall slope of the V φ -R curve making it flat, decreasing, and rising. We also note less ridges and less oscillations for the rising curve model (1hr) due to the slower winding rate, as already noticed in Fig. 4. The ridges in the third column appear more conspicuous and elongated in Model 1hd since the V R pattern, which follows lines of negative slope in this diagram, is favoured by the decreasing circular velocity curve. In fact, the V R pattern in the third column follows lines of constant angular momentum (in accordance with what is seen in the fourth column). On the contrary, the slope of the ridges in the V φ -R (second column) changes with the slope of the circular velocity curve combined with the sawtooth pattern. As for the wave of V R with L Z , we note that the wavelength varies between models with different n.
To examine the azimuthal dependency of the ridges, in Fig. 8 we show the radial velocity as a function of L Z and φ for different models with different circular velocity curves (columns) and times (rows). As already explained, for the kinematic spiral arms, the regions of negative (blue) V R correspond to the arms. The black lines show the locations of the arms in this projection. They were obtained by computing the minima of the V R wave following Eq. 7, where in practice one needs to plot L Z as a function of φ mod π. In the first row (t = 0), there is the perfect quadrupolar alternation of negative and positive V R regions for all models as imposed in the initial conditions. As time goes on, more wraps of the spiral are seen in the shape of inclined bands. As already seen, the number of wraps (bands) depends both on the time and the slope of the circular velocity curve. In the same manner, the (negative) slope of the blue and red pattern increases as a function of time and decreases with L Z at a fixed time.
As a final analysis with these simple models, we simulated two consecutive impacts which would reflect a more general situation during a satellite interaction with a disc. To do this, our A&A proofs: manuscript no. aanda initial conditions are those of Model 1hf after 2 Gyr of the first impact. At that time the particles receive a new velocity kick following Eq. 4 with ∆V 2 = 15 km s −1 , thus higher than the previous one (which was of ∆V 1 = 10 km s −1 ). The initial conditions and the configuration at 0.5 and 1 Gyr of evolution after the second impact are shown in Fig. 9. A series of largely wrapped spirals remaining from the first impact are observed at the initial stage to which the new pattern of Eq. 4 is superposed. At 0.5 Gyr a combination of the two patterns is visible, although the later one dominates visually due to its larger wavelength and amplitude. The inspection of the azimuthal sector at φ = π (Fig. 10) shows that the superposition of the two patterns is visible in the ridges and in the angular momentum space as two superposed waves with different wavelengths.
The inspection of the figures in this section indicate that the wavelength of the oscillations in various projections give information on both the time after perturbation and the circular velocity curve. We even see that signatures of multiple impacts may persist as superposed waves with different wavelengths. We therefore explore this further in Sect. 4, using Fourier analysis.

Frequency analysis
In this section we perform Fourier transforms (FT) of the signatures of the kinematic spiral arms in phase space to estimate the time of impact (i.e. the start of the phase mixing process) using the fact that they increase their frequency with time.
As already noticed, at a given time the frequency decreases with L Z , making the frequency of this signal ill-defined. Motivated by the shape of the curves of minimum V R in the L Z -φ space (Eq. 7) we can transform the L Z coordinate to make it lin-early dependent on φ: In this way, the curves of minimum V R are straight and equally spaced in the L n−1 n+1 Z -φ space. For the particular case of n = 0 in a flat circular velocity curve, this means simply using L −1 Z : As an example, in Fig. 11 we show the different flavours of Model 1 at time t = 1 Gyr in the transformed space L n−1 n+1 Z -φ. These depict straight and equally spaced bands. The black lines correspond to the analytical description given in Eq. 9. If a certain azimuthal coverage is available in the data, a fit of Eq. 9 to the blue bands gives, for known V 0 , n and R 0 , a direct determination of φ 0 and t from the slope and intercept. In practice one could also do a 2D FT for this fit. However, since our real data (see Sect. 6) do not cover a wide range of φ we explore how to gain (partial) information from data with limited azimuthal coverage.
From Fig. 11 and from Eq. 9 we can see that in a certain azimuthal location (φ) the valleys of minimum V R are now equispaced a distance (i.e. the wavelength): Therefore, with this equation we can try to recover the time of impact t in our simple models by measuring the frequency of the oscillations in V R with a FT in this transformed coordinate. We start by an example with Model 1hf in a given azimuth that we take as φ = π corresponding to the one in Fig. 6, and first choose the time t = 2 Gyr. The signal (V R oscillation) and FT (amplitude against frequency) are shown in the top and bottom panels of Fig. 12, respectively. We use V R /L Z to obtain a signal with constant amplitude with L Z (countering the imposed initial conditions, Eq. 4), although this does not affect the FT. For comparison, we do the FT of the L Z -V R curve (blue colours and bottom horizontal axis) and the case of the transformed coordinate L n−1 n+1 Z -V R , that for this case of n = 0 is L −1 Z -V R (orange curves and top horizontal axis). In the FT in the L Z coordinate we see multiple peaks over a wide range of frequencies. This is due to the fact that the frequency is not constant along the L Z axis. When we use L −1 Z , however, the signal shows a uniform frequency and the FT reveals a clear peak. Figure 13 (left) shows the amplitude of the FT now for all times for the same Model 1hf (colour map). A clear diagonal band is observed indicating an increase in the frequency with time following the expected ∝ t from the inverse of Eq. 11. The analytical formula is superposed in white and we see a perfect match between the predicted and the computed frequency, meaning that the FT analysis allows to recover the impact time for this simple model. In the two following right panels in Fig. 13 we show the computed frequency for the models with different slope of the circular velocity curve n. As expected, models with circular velocity curves with positive (negative) slope have a slower (faster) rhythm of phase mixing. We note that to do the FT in the transformed space L n−1 n+1 Z , we have assumed that we know n. Alternatively, one could search for the n that produces a clearer FT peak.  Finally, in the right panel of Fig. 13 we show the FT amplitude of Model 1hf with two superposed perturbations offset by 2 Gyr (see Fig. 10). For this idealised model, two lines corresponding to the signal of each individual perturbation are identified at all times, from which we can recover both impact times using Eq. 11.
In this section we have built simple models in which the dynamics can be described analytically. We now explore the more realistic N-body model using the simulation of Laporte et al.

N-body model
Our last model is the simulation L2 by Laporte et al. (2018). Briefly, this idealised live N-body simulation followed the interaction of a Sgr-like object with a MW-like host from the time of virial radius crossing to the present-day. The simulation was run with the Gadget-3 code (Springel 2005). The host dark halo and bulge were represented by Hernquist (1990) profiles with masses M h = 10 12 M and M b = 10 10 M , and scale radius of a = 53 kpc and a = 0.7 kpc, respectively. The disc followed an exponential profile in the radial direction and a sech 2 profile in the vertical one with scale radius R d = 3.5 kpc and scale height z d = 0.53 kpc, and a total stellar mass of M d = 6 × 10 10 M . The model was generated with GalIC (Yurin & Springel 2014) which was modified to include adiabatic contraction following the (Blumenthal et al. 1986) implementation. The Sgr-like progenitor was modelled as a Hernquist sphere of 8 × 10 10 M and a = 8 kpc, made to match a NFW halo of virial mass 6 × 10 10 M with concentration c = 28. This corresponds to a concentration that is about twice as concentrated as the mean of the Gao et al. (2008) mass-concentration relation which has about 0.3 dex scatter. For more details we refer the reader to (Laporte et al. 2018).
The pericentre passages of Sagittarius occur approximately at 2.3, 4.5, 5.6, 6.2 and 6.4 Gyr after the start of the simulation (the last pericentre roughly corresponds to the present time). We present the analysis of this model in four different time intervals starting always at a time previous to a pericentre: 1) the first peri-Article number, page 9 of 19 A&A proofs: manuscript no. aanda Z . The diagonal lines in all models show the increase in frequency with time, which is faster for the decreasing circular velocity curve mode (1hd). In the ideal model with two perturbations, two lines of the frequency corresponding to the two different impacts are detected at all times. centre (Model L18p1), 2) the second one (Model L18p2), 3) the third one (Model L18p3), and 4) the final sequence of snapshots when the last two pericentres (fourth and fifth) occur (Model L18p45). There are notable differences between these four time intervals. First, we examine the disc structure seen face-on and the velocity moments in Figs. 14, 15, 16, and 17, corresponding to the different time ranges.
Immediately after the first pericentre (panel at 2.4 Gyr in Fig. 14) there is a situation very similar to our simpler models: a quadrupole pattern in velocity (second and third columns) closely following the behaviour described by Eq. 4 (but already slightly coiled at 2.4 Gyr). Indeed, as shown in Fig. 11 of Grion Filho et al. (2021) where the same model is explored, there is a strong m = 2 mode in the Fourier decomposition of the radial and azimuthal accelerations at 2.2 Gyr (together with an m = 1 mode) for particles at radius of 10 kpc. With time, the velocity patterns wrap and a bi-symmetric spiral structure in density appears (first column). As before, the spiral arms correspond to regions of negative radial velocity and null azimuthal velocity with respect to the azimuthal mean. A dipole in V φ appears in the very central kpc which might be due to small centring uncertainties.  For the next time interval, in the first row of Fig. 15, we see that at 4.4 Gyr (that is 2 Gyr after the first pericentre) the spiral wraps from the first pericentre are still visible both in density and in velocities. After the second pericentre (panel at 4.6 Gyr), we again see a dominant m = 2 mode in the velocities that once more starts winding and forming a two-armed spiral structure. The impact of this second pericentre is slightly stronger than the first one, reflected in larger kick velocities (for example comparing panels at 2.4 and 4.6 Gyr). The signatures of the small scale ripples from the previous pericentre are hard to distinguish soon after the next impact (hardly visible by 0.3 Gyr after the second impact; middle row of Fig. 15). However, the signal of the previous impact can be seen better in other phase space projections (see below).
After the third pericentre (Fig. 16), a similar generation of a two-armed pattern is observed in density and in the velocity moments. In this case, though, the initial velocity perturbations follow Eq. 4 (quadrupole signature) only in the inner disc but not in the outer parts. This is due to the perturber crossing at a much closer distance (in fact it crosses the disc) and thus the condition of a distant impact is only roughly fulfilled in the inner parts. The resulting inner spiral, especially in velocities, is not as symmetric as before but in general evolves as in previous pericentres. The spiral structure in this case is even stronger than previously (both in density and velocities), as expected after a closer pericentre. Again some ripples stemming from the previous (second) impact are visible at the first times after the third pericentre. The outer parts of the disc follow slightly different patterns in velocity. The time span in Fig. 17 includes the last two pericentres of this model. These occur at separations shorter than the dynamical times of the stars in the disc outskirts (Grion Filho et al. 2021), producing large impact on the disc dynamics. The density structure of these final stages of the simulation is more complex, with two spiral arms but also rings and a central bar with ansae (appearing at about 6.1 Gyr). The velocity patterns associated to the spiral arms are not so clear, perhaps masked by the signal of the bar. The inner quadruple of the bar is obvious and the quadrupole with a phase shift of π/2 in some snapshots (e.g. 6.3 Gyr) could be due to the surpassing of corotation (Mühlbauer & Dehnen 2003) or the start of two new tidal arms given that these structures appear somehow sheared in the following snapshot. Now we examine the dynamics corresponding to a certain azimuthal angle in the disc (φ = π ± 0.2) examining the four different time ranges in Fig. 18, 19, 20 and 21, as we did for the simple models (Fig. 6, 7 and 10). For this model, we could not use Eq. 6 and thus in these figures we visually located the position of the overdensities in the histogram (first column) that correspond to the arm locations (or the bar in some cases) and marked them with vertical lines. For the vertical lines in the L Z -V R projections (fourth column) again we could not use Eq. 7 but we estimated the location of the minima by doing L Z = RV c as in the analytical model but with V c obtained from the N-body model, in particular, at a time previous to the perturber.
Starting with the first pericentre (L18p1), in Fig. 18 we observe ridges and an V R wave with a very similar behaviour to  Models 1. This is not surprising since we already saw similarity in the X-Y projection. Ridges form in a sawtooth shape (second column), with the ascending parts corresponding to the spiral arms in density (marked with vertical lines) and negative radial velocity, and positive radial velocities in the interarm regions and the descending ridges. We again see alternation of positive and negative radial velocities in the third column 2 . By 3 Gyr, three spiral arm wraps are seen in the histogram of radius (first column) with three associated ridges in the middle columns. A wave of V R in L Z forms with a progressively shorter wavelength with time similar to Models 1. The vertical lines in this projection correspond almost perfectly to the minima, as in our simple models.
The situation after the second pericentre (Fig. 19) is reminiscent of our simple model with two consecutive pericentres ( Fig. 10) but with the added complexity in the L18 model. Eventhough we saw that the small spiral ripples in the X-Y projection that remained from the first pericentre seemed to disappear in a relatively short time, here we keep seeing the corresponding ridges: at 4.6 Gyr just after the new pericentre occurs, the velocity V R seems overall dominated by negative values (blue) but we T. Antoja et al.: Tidally induced spiral arm wraps encoded in phase space  still see some small scale pattern of blue-red alternate (as evidenced also by the wave in L Z -V R in the fourth column). By 4.8 Gyr, for example, one clear bump in the histogram of R is seen corresponding to the new spiral wrap at this azimuth (clearly seen in Fig. 15) but we observe ridges with substructure in the second column and multiple blue and red bands in the third column. In these panels some of the thinner ridges do not follow the sawtooth shape anymore. This is because the two patterns of ridges (from the current and previous pericentre) superpose.
Overall the L Z -V R wave in Fig. 19 follows the same behaviour as in the earlier pericentre but with additional higher frequency undulations. The effects of the first pericentre are thus still visible after 3 Gyr after the first impact (at 5.2 Gyr, bottom panels). Although we suspect that this substructure also exist in configuration space, it is much harder to see due to the statistical noise. As in other work in different context, the velocity substructure seems to be preserved much better (Helmi & White 1999). Finally, the vertical lines in the fourth column coincide only approximately with the minima of the wave. This could be because in this case it is more difficult to locate the spiral arms in the first column due to the interference between the two perturbations and because the use of the V c from early times is not accurate anymore. The dynamics of the third pericentre (Fig. 20) also shows waves of different wavelengths from the different pericentres, similarly to L18p2.
In the final 2 pericentres (Fig. 21) there are clear bumps in the R histograms (associated to the spiral structure) that correlate with the thicker up-going ridges in the second column and with the negative radial velocities bands in the third column (and valleys in the fourth column). However, the correspondence between vertical lines and minima is good only in some cases possibly for the same reasons stated above. The spiral structure and associated dynamics seem more complex, probably due to the accumulation of effects of the multiple (and progressively stronger) impacts and the formation of the bar and rings, as already seen in Fig. 17. Figure 22 compares the azimuthal variation of the radial velocity wave in angular momentum at different times of the L18 model. Columns are the four different time spans specified at the beginning of the section. Overall the maximum velocities are larger for the last pericentre, as expected. The first pericentre is dominated by two approximately flat lines of alternate V R initially 3 that progressively become stepper with time and for smaller L Z , very similar to Models 1. The evolution after the second and third pericentres (second and third columns) globally follows a similar behaviour with hints of the superposition of the previous wave which create substructure inside each of the negative and positive regions. At larger angular momentum we see structure not following the diagonal patterns. This could be due to larger effects in these outer disc parts due to the new encounter (not fulfilling exactly the distant tide impulsive approximation). In addition or alternatively, it could be because the signatures of the first pericentre are less phase mixed due to the longer orbital timescales, in practise superposing waves of similar wavelengths. In the last time range, the kinematics at small angular momentum seems affected by the bar, as already mentioned. At larger angular momentum, there are more complex colour patterns than the inner diagonal stripes seen at other times and L Z . Now we apply the FT to the signal in this simulation pursuing the measurement of the increase in frequency with time after the different perturbations. We arbitrarily focus on the region around φ = π as in previous examples. Figure 23 shows the amplitude of the FT for all times. In the top panel we show the FT done in L Z space and on the bottom we use L n−1 n+1 Z , motivated by the discussion of Sect. 4 and the analytical formula Eq. 9. We note that this approach assumes that the potential is of the power-law type with a well defined n (constant with R), while this is not exactly the case for the model L18 (see for example Grion Filho et al. 2021 for the circular velocity curve and frequencies). This is, however, a first exploration of the potential of this methodology and not a far off approximation to realistic Galactic potentials despite its simplicity. One could also do it by using the exact frequencies of the model, which is something we can try in the future. Here we explored different values of n between -0.1 and 0.1 and found that, globally, n ∼ −0.1 retrieve higher FT amplitudes. At certain times, in particular towards the end of the simulation, n ∼ 0.1 also returns large amplitudes. This could be due to a change of the orbital frequencies at later times. Here for simplicity we show as n = −0.1 (bottom panel in Fig. 23).
In Fig. 23 an increase in signal can be seen coinciding with the pericentres (marked with white vertical lines) in both panels, but the signal seems clearer when the FT is done in L Z . This could be a sign that we are not using the proper transformation to L Z as explained above. This makes the process of age dating the perturbations from the analytical formula difficult. The poor signal might also be due to the resolution of the model. In any case, it is still possible to see diagonal bands of frequency increasing 3 We see again some short wavelength pattern that remains from the relaxation of the initial conditions as explained in footnote 2. linearly with time and starting at the first two pericentres. The FT amplitude also increases at times close to the last pericentres but we do not distinguish individual bands probably because of our time resolution and because the pericentres occur now very close in time and in the regime where bar effects start appearing.
An interesting question is whether we can detect multiple perturbations with different onset times simultaneously. While some signal could be seen directly in Fig. 18-21 and 22 (in some cases well after successive impacts), this is not so clear in Fig. 23 where the start of a new diagonal line seems to wash out the signs of the previous one. At this point it is not clear whether this is a real effect due to the later perturbations being stronger, or it is a detection problem since we could not find an optimal transformation of L Z and/or binning problem. In the future one could explore doing the FT in 2D in the full L Z -φ space, which most likely increases the signal, as pointed out in Sect. 4.

Data
The ridges in the R-V φ projection have been studied extensively but for an easier comparison with previous models we show again some figures with real MW data using radial velocities and astrometry from DR3 (Gaia Collaboration et al. 2022) and distances from StarHorse (Anders et al. 2022). We also applied the following quality cuts: i) astrometric quality selections, RUWE < 1.4 & parallax_over_error > 5, and ii) selection of non-spurious solutions (Rybizki et al. 2022), fidelity_v2 > 0.5. We also removed stars with rv_template_teff ≥ 8500 K since these stars present a residual bias of a few km s −1 in radial_velocity even after some proposed corrections (Blomme et al. 2022). In addition, we applied a correction to the radial_velocity to stars with grvs_mag ≥ 11 and rv_template_teff < 8500 (Katz et al. 2022) for which we use radial_velocity − 0.02755 × grvs_mag 2 + 0.55863 × grvs_mag − 2.81129.
For Fig. 24 we selected a wedge of φ ± 0.2 rad around the Sun-Galactic Centre line. The density structure in Fig. 24 (which is now in logarithmic scale) is the result of a complex convolution of the sample magnitude distribution and the different types of stars (most of the dwarf stars being located close to the Sun) with the true underlying density profile of the Galactic disc. In the second and third columns we see multiple thin and close ridges, resembling the later stages of a phase mixing process better than the initial ones. The existence of four minima in V R as a function of L Z (third panel) would indicate four wraps of spiral arms in the observed radial range if these were all due to a unique satellite impact but no other perturbation. We see a fifth minima at low angular momentum but at this range of L Z selection effects are more dominant (favouring eccentric orbits) and this position is closer to the region of the bar and could thus more likely be unrelated to the phenomena that we study here.
The ridges seen in the real data do not show such a clear sawtooth shape in Fig. 24 as in some models presented in Sects. 2-5. However, this might depend on the particular plotting technique or sample used since some zigzagging seems present in Fig. 12 in Hunt et al. (2019) andFig. 16 in Gaia Collaboration et al. (2021). We also note that the sawtooth shape is not as strong T. Antoja et al.: Tidally induced spiral arm wraps encoded in phase space  in the simple model with decreasing circular velocity curve (as perhaps in the MW) nor in the more realistic L18 model. The sawtooth shape also disappears when multiple perturbations coexist, which is very likely the case for our Galaxy. Additionally, the zigzag is present also in models that do not have tidally induced arms, such as those by Khanna et al. (2019) and Hunt et al. (2018). Therefore, this shape can not be used to disentangle between models.
A wave in the L Z -V R space was already detected in Friske & Schönrich (2019) who measured two wavelengths of 1350 and 285 km s −1 kpc and attributed the larger one to the signatures of the bar. The new, larger sample used here shows a similar result (right panel in Fig. 24 and Fig. 25) but much better defined oscillations and larger radial (L Z ) and azimuthal span. Following up on previous sections, we also compute the FT for the data. One could try to do a fit of Eq. 9 directly in the L Z -φ space. However, the range of φ is small and we do not see a regular slope of the coloured bands in Fig. 25. This is already a warning of the higher complexity of the data compared to our models. Yet we do the FT to the L Z -V R curve at φ = φ , and attempt to do a fit using Eq. 11 for a back-of-the-envelope calculation. To do this, we need to assume an n (slope of the circular velocity curve). We considered different circular velocity curves of the MW in the literature, namely from McMillan (2017), Sofue (2020), Eilers et al. (2019) and MWpotential2014 from Bovy (2015), and we find the n that best fits each curve in the range where we detect our signal (5.5-12 kpc). We obtain n =0, -0.03, -0.06, and -0.1, respectively. We note, however, that only the curves by Eilers et al. (2019) and MWpotential2014 are well fitted by a single power law model. By contrast, McMillan (2017), for example, shows raising and decreasing trends with the limit at R ≈ 8 kpc.
In Fig. 26 we show two curves corresponding to the V R signal (top) and FT amplitudes (   (i.e. L −1 Z ), as an example from all the above circular velocity curves. We find that L Z gives slightly more defined frequency peaks than L −1 Z and again this could be due to a bad approximation of the real potential to the power law type. In this case we also use V R /L Z but it does not affect the results of the FT.
In the L Z coordinate (blue) we find two defined peaks at a small and a large frequency. We note though that the peak at smaller frequency is at the limit of the first frequency point allowed by the data, and thus only an upper limit can be measured. Since our sampling of the frequencies is sparse, we only give a range of possible location of the peaks. We measured this range with the half-distance between the peak the consecutive frequency points. These error ranges are marked as horizontal blue error bars on the bottom panel of Fig. 26. The frequencies that we find correspond to wavelengths of > 1100 and in the range of 290-360 km s −1 kpc, thus similar to Friske & Schönrich (2019). For L −1 Z and in the case of the larger frequency zone, a broad flat peak is seen at approximately 10 000 (km s −1 kpc) −1 . Again, for the smaller frequency we can only measure an upper limit.
The limiting frequencies obtained in each case are given in Table 2. In all cases, we computed the times of impact using Eq. 11 and assuming V c = 239.26 km s −1 and R 0 = 8.277 kpc. As expected, the obtained times are sensitive to the adopted value of n, with larger n yielding larger times. Considering all cases, our analysis gives times of around < 0.6 Gyr from the low frequency peak and of 0.8-2.1 Gyr from the larger frequency one. We discuss these estimated times in Sect. 7.2 .

Spiral arms, ridges, waves, and frequencies
The response associated to the density wave theory in the tightwinding approximation (TWA) predicts that density and average radial velocity are anti-correlated inside corotation radius (i.e. the spiral arms coincide with negative radial velocities) but are positively correlated outside corotation (Lin et al. 1969;Mayor 1970). The TWA arms also coincide with regions of null azimuthal velocity with respect to the average at that radius (as expected from the relation between radial and azimuthal velocities in disc orbits). In Antoja et al. (2016) (see also Siebert et al. 2012 andRoca-Fàbrega et al. 2014) we found that this was indeed approximately the case in test particle simulations with potentials including TWA arms with a constant pattern speed. However, the self-consistent isolated N-body models examined in that study with strong spiral arms and with transient arms showed more complexity. More recently, in Eilers et al. (2020) the kinematic perturbation associated to logarithmic spiral arms was calculated through approximating orbits in a potential including the spiral perturbation. The resulting peaks of highest density, i.e. the spiral arms, coincided with locations of zero average radial velocities. More examples of the analysis of the spiral arms' kinematics can be seen in Grand et al. (2015) and Pettitt et al. (2020) who also focused on the density wave theory, transient spirals, and spirals coexisting with bars. The spiral structure of our models induced by an external perturber corresponds to phase-space regions with negative average radial velocities and zero azimuthal velocity with respect to the average. This is consistent with the theoretical considerations of Kalnajs (1973) about kinematic waves and was observed in early N-body modelling of a galaxy perturbed by a point mass (Sundelius et al. 1987) and later, for example, in Oh et al. (2008).
Here we put this in the context of spiral arms that could be tidally induced by the successive encounters of the MW and Sagittarius. This clear kinematic signature can potentially help to distinguish different types of spiral arms described above.
Moreover, the tidally induced spiral arms correspond to diagonal ridges with negative radial velocity V R in the R-V φ projection, which does not necessarily correspond to regions of higher stellar density in this projection. Our spiral arms are also linked to the L Z -V R wave, corresponding to its valleys (negative V R ). Since the spiral arms in density space are much more difficult to detect in real data due to extinction and selection effects, our findings offer a potential method to infer the number of spiral arms and their location. This is different from Khoperskov et al. (2020), where the spiral arms were associated to regions of high density in the R g -φ plane (essentially regions of accumulation around certain angular momentum). Their spiral arms, thus, do not correspond to spiral arms necessarily seen in configuration space (see also Hunt et al. 2020), where we believe spiral arms must be defined. From our models we see that overdensities in certain projections of phase space do not correspond to overdensities in other projections, in particular not in configuration space.
In our models the number of ridges in the R-V φ plane, the number of alternating changes of sign of V R in this projection, and the number of oscillations in the L Z -V R wave (all these corresponding to spiral wraps) depend on the form of the planar potential and the time after the perturbation. For simple models with power law potentials, we give analytical descriptions of the number of wraps, positions of the ridges, etc, as a function of time t (after the perturbation), the initial phase φ 0 , and the parameters of the circular velocity curves (slope n and normalisation values V 0 and R 0 ). By measuring the frequency of the L Z -V R wave and assuming a reasonable model for the gravitational potential, one can obtain the time of onset of the perturbation. Ideally, one would try to fit the potential parameters and time of impact simultaneously. Here we also show that planar kinematic disturbances can contain signatures of multiple events as superimposed waves of different wavelengths. Our Fourier analysis allows the approximate detection of different impact times for these simple models, while there seems to be a certain detection limit for the analysed N-body simulation, in which there is added complexity due to other perturbations, a potential deviating from simple power laws, and low number of particles.
We note that there is a clear analogy between this and the estimation of the impact time using the phase spiral in Antoja et al. (2018) and also in using this phase-space structure to constrain the potential (Widmark et al. 2021), but our observable in this case relates to the planar potential and not the vertical one. As discussed in Laporte et al. (2019) and Hunt et al. (2021) for the phase spiral, multiple perturbations from the planar velocities could potentially be easier to identify in the outer parts of the Galaxy due to longer phase mixing times (Fig. 5). Developing a framework in which both the planar and vertical dynamics are modelled together is necessary and could provide a more robust fit to both disturbances at the same time. We now proceed to discussing the applicability of our idea to currently available Milky-Way data.

The case of the Milky Way
Our simple models with orbital integrations provide a basis to understand the dynamics of tidally induced arms and their manifestations in different phase-space projections. They have proven to be very helpful to understand what happens in more complex models such as that of L18, and potentially real disc galaxies such as the MW. One of the simplifications of our simple models is the condition of distant and impulsive impacts. A natural extension of our work could be to relax the distant-impact condition and explore different impact configurations in our toy models that could resemble the later interactions with a Sagittariuslike system, i.e. when the perturber may even cross the disc, and excite more complex structures such as asymmetric rings and spirals (Appleton & Struck-Marcell 1996). The other fundamental assumption in our simple modelling is the neglection of the effects of self-gravity of the disc and other perturbations, mostly the bar (and possibly other types of spiral arms). More investigation is needed in order to see how these additional aspects can hamper our ability to constrain the time of the perturbations and potential of the MW. For example, Darling & Widrow (2019) show that phase-mixing timescales are longer when one considers self-gravity. Similarly, Pettitt & Wadsley (2018) and Oh et al. (2015) show how the phase-mixing rhythms of tidally induced arms are altered by self-gravity or by the bar.
As both the existence of the bar and of the external perturber Sagittarius are well established, the Galactic disc must be suffering from sequels of both, the key questions being whether one of the two is the dominant mechanism, whether the two have superposing effects -and are thus detectable as separate waves -or whether they interact with each other. Hunt et al. (2019) studied the combined effects of a bar and spiral arms with test particles and backward integration of orbits, but they considered only quasi-stationary density wave arms and transient corotating arms. In a future study, our simple modelling could be used to examine the combination of bar and tidally induced spiral arms, complemented by a more detailed study of the later stages of the model L18 or other N-body models. In particular, one could explore how the FT analysis presented here could help to distinguish waves/signatures from different dynamical processes, such as the bar.
In the N-body simulation that we use here, the effects of the Galactic bar start to become important during the last pericentres, complicating the interpretation of the density and velocity structures due to the interaction. The kinematics associated to a barred potential has been well studied in the literature. For a recent example, we refer to Fig. D.1 of Bernet et al. (2022), which shows the R-V φ projection coloured by V R for test particle simulations with different bar pattern speeds at 30 deg with respect to the bar's long axis. We see that there are clear differences between the bar and the external perturbation model. As expected, there is a limited number of ridges in the bar models: the most clear ones are those associated to corotation, Outer Lindblad Resonance (OLR) and 1:1 resonance 4 . These ridges do not have the conspicuous sawtooth shape seen in our models, and some of them have changing slopes or breaks with R, and their V R signal is different depending on the resonance. This translates into a wave in the L Z -V R space that does not have the symmetries that we see in the waves of external perturbations in the simple models.
For the bar case, the V R wave in the L Z -φ plot (similar to Figs. 8 and 23) can be seen in Fig. A1 of Trick et al. (2021) and Fig. 12 of Chiba et al. (2021). There we see that the bar mostly creates a quadrupolar structure in V R that changes sign at the main resonances (see also Mühlbauer & Dehnen 2003). Trick et al. (2021) show that for long integration times, this structure has phase angle φ constant with L Z , thus different from the external perturbation case. Only for short integration times does one see certain slopes, but this could also depend on the resonance, thus different from the global pattern in our models. These aspects should be key to consider in the comparison with the Milky-Way data (see below). However, Chiba et al. (2021) show that this quadrupolar pattern turns into multiple bands and slopes when the bar slows down, which is a natural process in galactic dynamics. In this case, the simulations resemble more the data although additional effects might still be at play.
Finally, we note that simple barred models predict density gaps and bumps around some resonances such as the OLR (for example Fig. 9 in Trick et al. 2021, or Figs. 6 and7 in Melnik et al. 2021). While it is well described that the different bar resonances can create structures in density (Buta & Combes 1996), this has been hardly considered for the MW. These structures have mostly ring shape but we note that, locally, they might be confused with spiral arms. It should be explored further whether their kinematics makes them distinguishable from spiral arms.
In any case, awaiting a clearer picture of how different dynamical effects (bar, external perturbation, etc) superpose, our simple modelling indicates two times of impact of < 0.6 and 0.8-2.1 Gyr, respectively for the small and large frequencies detected in the data. These frequencies could in principle also be related to the effects of a static bar (Mühlbauer & Dehnen 2003;Chiba et al. 2021;Trick et al. 2021) in the small frequency case, as discussed in Friske & Schönrich (2019), and to a slowing bar according to Chiba et al. (2021) in the large one. From our simple modelling, the small-frequency signal (related to the short time) could alternatively be linked to a recent passage of Sagittarius, although this is expected to happen a bit earlier in time and to have conditions far from the distant-impact hypothesis. The large-frequency peak yields times compatible with the perturbation from the previous-to-last (or earlier) pericentre (e.g. de la Vega et al. 2015;Laporte et al. 2018;Law & Majewski 2010;Purcell et al. 2011). Signatures of previous pericentres would not be washed away if Sagittarius significantly lost mass between pericentres. Minchev et al. (2009) obtained a time of perturbation of 1.9 Gyr from the visual comparison of the separation between arches in the V R -V φ plane in the data and in their toy models. Their method is considerably different from ours, since we start with initial conditions that specifically mimic the ones following an impulsive and distant impact with a satellite. Furthermore, we derive an analytic formula for the oscillations in the (transformed) L Z -φ space, consider the azimuthal dimension of the phase-mixing, and relate the perturbation with the tidally induced spiral arms. In any case, one of our derived perturbation times is consistent with the one obtained by Minchev et al. (2009).
Our derived lower limit of the impact time related to the large-frequency signal (0.8 Gyr) is also comparable to the times inferred from the phase spiral in Antoja et al. (2018) (500-900 Myr) in the pure phase-mixing case (although it was estimated to be 400 Myr in Binney & Schönrich 2018), but this is expected to be larger in the self-gravitating case (Darling & Widrow 2019). The times are also compatible with direct age-dating of the phase spiral (Tian et al. 2018;Laporte et al. 2019;Bland-Hawthorn et al. 2019), which seems to be present in young samples (0.5-1 Gyr).
We note that the amplitude of the V R wave in the last pericentres of the L18 model is too large compared to the data ( 20 vs 10 km s −1 ). This could simply be due to the lack of exploration of different setups for Sgr's halo and disc structure in the current models. It could also be connected to the existent discussion around the mass of Sagittarius at each pericentre and whether these impacts are enough to create the phase spiral (see discussion in García-Conde et al. 2022). This should also be addressed in future investigations.

Summary and conclusion
We have presented a series of models of tidally induced spiral arms and shown how these relate to the ridges and waves seen in certain projections of phase space, similar to those observed in the Gaia data for the MW. Our main results and conclusions are: 1. An impulsive distant tidal impact generates an initial quadrupole in velocities that leads to a two-armed spiral structure that winds up with time at an approximate rate of ω − 1/2κ. The arms progressively have smaller pitch angles and there is an increasing number of spiral wraps in a particular azimuth in the disc with time (Sect. 2.2). 2. The location of the maximum density (locus of the arms) coincides with regions where the average galactocentric radial velocity is negative, and the azimuthal velocity with respect to the average at each radius is zero (Sect. 2.2). This is because the arms coincide with orbital phases in-between apocentre and pericentre. 3. The winding up of the spiral arms and the associated velocity pattern produce a global wave in V R that depends on azimuth, manifesting as inclined stripes in the L Z -φ plane, and as oscillations of V R as a function of L Z in a given region of the disc (Sect. 3), as observed also in the Gaia data. The wavelength of this oscillation decreases with time (and with L Z at a fixed time). This wave is seen as a pattern of stripes with alternating signs of V R in the R-V φ projection. Each stripe with negative V R corresponds to a spiral arm wrap. 4. In the R-V φ plane, a sawtooth shape appears, with increasing number of oscillations with time (Sect. 3). Each oscillation is associated to a spiral wrap, but the stars that belong to the arms do not lie on the down-going ridges but on the regions with negative radial velocity in this projection. The sawtooth shape is diluted in more complex models such as N-body and when multiple mechanisms are at play (Sect. 5). 5. The relation between density and average velocities (and associated ridges and waves in phase space) is different for other models of spiral arms, offering a way to distinguish between different dynamical models of arms (Sect. 7.1). 6. With the realistic simulation of the interaction between Sagittarius and the MW, we see that the first three pericentres of Sagittarius seem close to the impulsive and distant tide approximation in the central and intermediate regions of the disc, and thus, the dynamics follows the items above (Sect. .5). 7. The rhythm at which new spiral arms wraps appear with time, and thus new ridges in the R-V φ projection and oscillations of V R with L Z , depends on the orbital frequencies of the disc through the radial dependence of ω − 1 2 κ. The number of all these phase-space features and their separation offers a means to infer at the same time the potential (essentially the circular velocity curve parameters) and the starting time and phase of the perturbation. We give analytical formulae for simple power law models that at a theoretical level (in pure phase mixing events) successfully recover the impact times and even disentangle impacts from different successive pericentres (Sect. 4). 8. For the MW we find a superposition of two different wavelengths of > 1100 and 290-360 km s −1 kpc, similar to the values of Friske & Schönrich (2019) who detected it for the first time (Sect. 6). 9. Using typical values of the circular velocity curve of the MW, we obtain times of < 0.6 Gyr and 0.8-2.1 Gyr, associated to these two wavelengths (Sect. 6). The first one is slightly smaller than the last Sagittarius pericentre and thus the long wavelength signal could be related to other mechanisms. The range of earlier times is consistent with estimated previous and previous-to-last pericentre of this galactic satellite. These would correspond to four spiral wraps in the range of the data coverage. 10. The validity of this calculation is subject to the assumptions that the impact with Sagittarius was impulsive and distant, and that there is no intervention of other aspects such as selfgravity or the bar, which, on the other hand, has been shown to reproduce several of the observed kinematic features elsewhere as well (Sect. 7.2).
Here we provided various models of the dynamics of tidally induced arms and a first application to the Milky Way case. However, the complex dynamics of the MW disc certainly requires more investigation.