# Measuring the vertical response of the Galactic disc to an infalling
satellite

Eloisa Poggio${}^{1,2},$ Chervin F. P. Laporte3,4,2, Kathryn V. Johnston5,6,
Elena D’Onghia7, Ronald Drimmel1 and Douglas Grion Filho5
1Osservatorio Astrofisico di Torino, Istituto Nazionale di Astrofisica (INAF),
I-10025 Pino Torinese, Italy
2Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, France
3 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study (UTIAS),
The University of Tokyo, Chiba 277-8583, Japan
4 Department of Physics & Astronomy, University of Victoria, 3800 Finnerty
Road, Victoria BC, V8P 5C2 Canada
5Department of Astronomy, Columbia University, 550 West 120th Street, New
York, NY 10027, USA
6Center for Computational Astrophysics, Flatiron Institute, 162 5th Av., New
York City, NY 10010, USA
7University of Wisconsin, Madison, Astronomy Department, 475 N Charter Str.,
Wisconsin, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Using N-body simulations of the Milky Way interacting with a satellite similar
to the Sagittarius dwarf galaxy, we quantitatively analyse the vertical
response of the Galactic disc to the satellite’s repeated impacts. We
approximate the vertical distortion of the Galactic disc as the sum of the
first three Fourier azimuthal terms $m=0,1$ and $2$, and observe their
evolution in different dynamical regimes of interaction. After the first
interaction, the $m=0$ term manifests itself as outgoing ring-like vertical
distortions. The $m=1$ term (S-shape warp) is prograde when the impacts of the
satellite are more frequent, or in general close to an interaction, whereas it
is slowly retrograde in the most quiescent phases. The $m=2$ term is typically
prograde, and close to an interaction it couples with the $m=1$ term. Finally,
we find that the vertical response of the disc can be recovered in an unbiased
way using the instantaneous positions and velocities of stars in a limited
volume of the Galactic disc, analogous to real data, and that the measured
vertical pattern speeds have a constraining power in the context of a Milky
Way-satellite interaction.

###### keywords:

Galaxy: disc – Galaxy: evolution – Galaxy: kinematics and dynamics – Galaxy:
structure - galaxies: Local Group

††pubyear: 2015††pagerange: Measuring the vertical response of the Galactic
disc to an infalling satellite–C

##  1 Introduction

The Milky Way’s vertical structure may hold important clues of its past
formation. As soon as 21-cm surveys made it possible to study the distribution
of neutral hydrogen gas around edge-on galaxies, it became apparent that many
galaxies, including our own, are warped (Burke, 1957; Kerr, 1957). With the
advent of large photometric surveys, such as the 2MASS and Sloan Digital Sky
Survey, it became evident that the disc exhibited more complex behaviour in
its vertical structure. In the solar neighbourhood, Widrow et al. (2012)
uncovered signs of vertical waves in the form of North/South asymmetries in
number density and bulk velocity of stars. The kinetic signatures of these
waves were confirmed in Williams et al. (2013); Carlin et al. (2013);
Schönrich & Dehnen (2018). Moreover, a growing number of observational studies
suggested these waves could be part of much larger scale ones, giving rise to
a warped and rippled disc (Xu et al., 2015; Price-Whelan et al., 2015; Li et
al., 2017).

A number of mechanisms have been proposed to explain the existence/persistence
of warps in galaxies. These have ranged from special configurations (and
assumptions) such as the presence of flattened potential around discs (Sparke
& Casertano, 1988; Dubinski & Kuijken, 1995), torques exerted by rotating
halos (Nelson & Tremaine, 1995), misaligned angular momenta between halo and
disc (Debattista & Sellwood, 1999) to misaligned gas infall (López-Corredoira
et al., 2002) and the influence of external perturbers (Weinberg, 1998), in
particular through the excitation of global dark matter wakes which can
penetrate deep into the Galactic midplane as well as tidal interactions. In
particular satellite interactions are a prime mechanism to producing not only
warps (e.g. Gómez et al., 2017; Semczuk et al., 2020), but also the more
complex structure such as the generation of bending waves (e.g. Gómez et al.,
2013; Widrow et al., 2014; D’Onghia et al., 2016), rings (Younger et al.,
2008; Kazantzidis et al., 2009) and kinematic features as tidal tails (Toomre
& Toomre, 1972; D’Onghia et al., 2010; Laporte et al., 2019a) which are seen
in the Milky Way. In fact, tidal encounters are central to the lifetime of
galaxies in a cosmological context, as every central galaxy is surrounded by a
subhalo population which can in some cases (depending on the mass and distance
reached by a subhalo) shape the evolution of Galactic discs (e.g. Font et al.,
2001; Gauthier et al., 2006; Aumer & White, 2013, Giammaria et al.,
submitted). Given the characteristic functional form of the subhalo mass
function (Springel et al., 2008), the less numerous but much more dynamically
impactful subhalos are the more massive ones.

In the Milky Way, potential suspects have long been the Magellanic Clouds
(Weinberg & Blitz, 2006) and the Sagittarius dwarf (Ibata & Razoumov, 1998;
Bailin, 2003) galaxy (though with mixed success until more recently). While
the Magellanic Clouds are more massive, they are also more distant, and new
proper motion measurements (Kallivayalil et al., 2013) indicate that they are
on a first infall, implying only a very recent interaction. However, first
infall N-body models of the Large Magellanic Cloud (Laporte et al., 2018a)
show that they are still able to warp the Galactic disc similarly to the HI
observations (Levine et al., 2006) through the excitation of a wake which may
be traced using the stellar halo (e.g. Garavito-Camargo et al., 2019). On the
other hand, Sagittarius (hereafter, Sgr) is currently undergoing full
disruption, and has been perturbing the Galaxy for a more extended period of
time through numerous pericentric passages. Moreover, recent estimates of its
total stellar progenitor mass place it at higher infall halo masses (e.g.
Gibbons et al., 2017) which significantly raises its potential to shape the
Galaxy. In particular, it has been proposed that Sgr could be an architect for
spiral structure in the Milky Way (Purcell et al., 2011) and the formation of
bending waves seen in the solar neighbourhood and at large distances (Gómez et
al., 2013) at least in a qualitative sense. Recent N-body simulations of Sgr-
like interactions with the Milky Way (Laporte et al., 2018b) show that it is
possible to simultaneously reproduce the amplitude of the density fluctuations
in the solar neighbourhood and the extent of known outer disc structures
(Newberg et al., 2002; Rocha-Pinto et al., 2004).

Thanks to the second _Gaia_ Data Release (hereafter, _Gaia_ DR2, Gaia
Collaboration et al., 2018a), the complexity in kinematics and richness in
terms of substructure of the Galactic disc has been mapped to greater detail
and unprecedented spatial coverage around the Sun (Gaia Collaboration et al.,
2018b). Vertical disturbances, streaming motions, wave patterns and arches in
the velocity space have been revealed on a large scale (Ramos et al., 2018;
Poggio et al., 2018; Romero-Gómez et al., 2019; Friske & Schönrich, 2019),
triggering new interest and questions into the physical mechanisms regulating
the dynamical evolution of the disc (Carrillo et al., 2019; López-Corredoira
et al., 2020; Poggio et al., 2020). Complementary to the earlier findings of
Widrow et al. (2012), when viewed in 2D, the vertical phase-space
distributions of stars in the solar neighbourhood revealed spiral snail-shell
features (Antoja et al., 2018), signaling that the disc is out of equilibrium
and currently undergoing phase-mixing, possibly caused by Sgr (Antoja et al.,
2018; Binney & Schönrich, 2018; Laporte et al., 2018b; Bland-Hawthorn et al.,
2019; Haines et al., 2019; Laporte et al., 2019b; Bland-Hawthorn & Tepper-
Garcia, 2020). Alternatively, phase-space spirals have been suggested as a
manifestation of internal mechanisms, such as the bar buckling instability
(Khoperskov et al., 2019). The sharper view of Gaia DR1 combined with the
SDSS’s long baseline to provide accurate proper motions, and later those of
_Gaia_ DR2, has also revealed remnants and confirmed earlier signs of ancient
massive accretion events, among which stands out the Gaia Enceladus-Sausage
galaxy (Meza et al., 2005; Navarro et al., 2011; Belokurov et al., 2018; Helmi
et al., 2018; Haywood et al., 2018; Recio-Blanco, 2018; Helmi, 2020), thought
to have occurred about 10 Gyr ago.

With the increasing complexity faced by the data and the number of competing
dynamical processes which could occur in concert, it is now imperative to
study and develop theoretical models to help with the interpretation of
observations. Ultimately, one would want to discriminate between the different
dynamical processes that might be at play (bar buckling, misaligned gas
accretion, halo rotation, satellite(s) interaction(s)…). Recently, the
vertical response of galactic discs to multiple external perturbers has been
revisited in some detail by Chequers et al. (2018), who compared the behaviour
of bending waves in the presence of a population of dark matter subhalos and
in isolation. However, the tidal effects are dominated by the most massive and
luminous substructures, and many subhalos never reach the central part of a
halo (Font et al., 2001). Thus, a complementary approach to Chequers et al.
(2018) is to study the response of a Galactic disc subject to a minor merger
over a lifespan of several orbital timescales on a polar orbit. This should
help get a better intuition of the effect of a particular kind of orbital
configuration shared by at least one massive known luminous satellite of the
MW, the currently disrupting Sgr dwarf.

In this work we aim to explore the reaction of a Milky Way-like disc to a
massive Sgr-like satellite as it spirals deeper into the galactic potential.
The forthcoming Gaia Data Release 3 and other future datasets will afford even
more detailed and global views of the Galaxy, motivating us to perform a
Fourier analysis of the disc vertical distortion, which may be feasible with
stellar data. In this paper, we aim to: (i) characterize the response of the
disc to the repeated impacts of the dwarf galaxy, (ii) study the connection
between the satellite’s orbit and the disc response; (iii) investigate whether
the response might be reliably measured from a limited data volume; (iv)
understand whether the vertical response of the disc has a constraining power
in the context of a Milky Way-Sgr-like interaction. In a companion work, we
look at the influences that drive this evolution (Filho et al., in prep.).

The paper is structured as follows. In Section 2, we describe the simulations
used in this study. In Section 3, we define and discuss different dynamical
regimes present in the MW-Sgr interaction process. Section 4 is dedicated to
the analysis of the disc vertical response to the repeated interactions with
the Sgr-like satellite. In Section 5, we compare the actual evolution of the
disc obtained in Section 4 to the results that would be obtained from a
volume-limited survey. Finally, we summarise and discuss our our results in
Section 6.

Figure 1:  _Upper panel_ : Vertical elevation of Sgr with respect to the
Galactic plane as a function of time, based on Sgr L2 model. _Middle panel_ :
Same as before, but showing the distance from the Galactic center. _Lower
panel_ : Regimes of dynamical interaction between Sgr and the Milky Way,
defined as described in Section 3. The grey horizontal line shows where
$T_{\rm{dyn}}/\Delta T_{\rm{impacts}}=1$, which roughly corresponds Regime 2.
In all three panels, the dashed vertical lines indicate Sgr’s plane-crossings
(i.e. where Z=0). Due to a different response of the disc (see text), we have
made a distinction between the plane-crossings in proximity (thick dashed
vertical lines) and distant from (thin dashed vertical lines) Sgr’s nearest
pericenter.  Figure 2: Face-on views of the Milky Way’s disc in the four
interaction regimes (see text: surface density (upper row), mean vertical
coordinate Z (middle row) and vertical velocity Vz (lower row) at different
times. At t=2.1 Gyr (first column) no impacts with Sgr have occured. At t=3.6
Gyr (second column), t=5.3 Gyr (third column), t=6.9 Gyr (fourth column), the
galactic disc has experienced, respectively, one, three and eight impacts.
Disc rotation is clockwise.

##  2 Simulations

The simulations analysed here are taken from the pre-Gaia DR2 runs performed
in Laporte et al. (2018b) (hereafter L18, ). These consist of a series of
controlled N-body experiments aimed at studying the response of a MW-like
system to Sgr and LMC-like perturbers. Despite its modest setup, this
simulation was used to establish that the large scale perturbations in the
outer disc such as the Monoceros Ring complex (Newberg et al., 2002; Slater et
al., 2014; Deason et al., 2018), TriAnd (Rocha-Pinto et al., 2004) and small
scale inner disc perturbations (Widrow et al., 2012; Antoja et al., 2018) can
be simultaneously understood as the response of the Galactic disc to
perturbations from a satellite like Sgr (Laporte et al., 2018b, 2019a). A
number of observaitonal studies had suggested a disc interpretation for outer
stellar halo-disc interface structures (Price-Whelan et al., 2015; Xu et al.,
2015; Li et al., 2017; de Boer et al., 2018; Sheffield et al., 2018) and were
later confirmed through chemistry (e.g. Bergemann et al., 2018; Hayes et al.,
2018). It was also the first N-body simulation used to interpret the Gaia DR2
RVS velocity fields and phase-space spiral in Gaia DR2 and study its formation
with a realistic model of a Sgr-like satellite interaction (Laporte et al.,
2019b). The simulations were also used to predict the formation of structures
like the Anticenter stream (Laporte et al., 2019a) as remnants of tidal tails
excited by satellite interaction and the expected modulation of the star
formation history of the Galaxy following the impacts which was first reported
in (Laporte et al., 2020) using the outer disc Monoceros complex and recently
in the solar neighbourhood (Ruiz-Lara et al., 2020).

In this paper, we take a closer look at the L2 run, which considered the
interaction between a MW-like host and a Sgr-like progenitor with
$M_{SGR200}=6\times 10^{10}\,\rm{M_{\odot}}$ infalling from the virial radius.
For more details on the simulations and their analysis, we refer the reader to
L18, but provide a brief summary here to guide the discussion. For all runs, a
fiducial host Milky Way-like model was adopted consisting of an exponential
disc in the radial direction and isothermal vertically with scale radius
$R_{d}=3.5\,\rm{kpc}$ and scale height $z_{d}=0.53\,\rm{kpc}$, a Hernquist
bulge with $M_{b}=10^{10}\,\rm{M_{\odot}}$ and scale length
$a_{b}=0.7\,\rm{kpc}$ surrounded by a dark matter halo of mass
$M_{200}=10^{12}\,\rm{M_{\odot}}$ approximated by a Hernquist (1990) profile
with scale length $a_{h}=53\,\rm{kpc}$ which is adiabatically contracted
according to the prescription of Blumenthal et al. (1986). This results in a
galaxy with a rotation curve consistent with McMillan (2011). The initial
conditions were realised with GalIC (Yurin & Springel, 2014) and the galaxy
was set with a Toomre $Q>1.5$ to avoid the formation of non-axisymmetric
features in the isolated case. This choice is arbitrary of course and only
intended for the sake of the numerical experiment at the time to understand
the response of the disc in the absence of internal perturbations from the
start.

Throughout the paper, we analyse the disc’s dynamics in the simulation from a
reference frame centered around the Galactic Center with respect to the
Galactic midplane. This is achieved using the same analysis methods outlined
in Gómez et al. (2013); Laporte et al. (2018b), which involve centering the
galaxy through the iterative shrinking sphere algorithm (Power et al., 2003)
and re-aligning the system with respect to the total angular momentum of the
inner disc. As a result, any re-orientation or precession of the disc w.r.t.
to its initial configuration resulting from interactions (e.g. Velazquez &
White, 1999; Gómez et al., 2017) will be lost. This dynamical response is not
studied in this paper.

##  3 Interaction regimes

The interaction process between Sgr and the Galactic disc can be divided into
different regimes, based on the the orbit of Sgr and the dynamical properties
of the disc.

At time t=0, the progenitor of Sgr is launched into orbit and crosses the
Galactic plane several times before merging fully with the Milky Way. As time
passes, plane-crossings become more frequent (Figure 1, top panel), with Sgr
passing through the Galactic plane (z=0 kpc) at 2.36, 4.43, 4.83, 5.58, 6.1,
6.31, 6.48, 6.73 and 6.95 Gyr (dashed vertical lines in Figure Figure 1). At
the plane-crossings, Sgr is located at 32, 30, 55, 11, 9, 24, 5, 2.5 and 1 kpc
from the Galactic center, respectively. As can be seen in Figure 1, most (but
not all) of the plane-crossings are close to the orbit pericenters.

As discussed in the following section, the influence of Sgr becomes
progressively more important as its orbit gets closer to the MW and the
impacts become more frequent. However, the response to Sgr is not the same
everywhere in the Galactic disc. We define the dynamical time $T_{\rm{dyn}}$
(i.e. Galactic year) as a function of Galactocentric radius R, i.e.
$T_{\rm{dyn}}\sim 2\pi R/V_{\phi}$, where $V_{\phi}$ is the mean azimuthal
velocity. For example, for stars rotating at $V_{\phi}\sim 230$ km/s, the
dynamical time would be $\sim$ 220 Myr at $R=8$ kpc, or $\sim$ 550 Myr at
$R=20$ kpc.

We define different Sgr-MW interaction regimes based on the ratio
$T_{\rm{dyn}}/\Delta T_{\rm{impacts}}$, where $T_{\rm{dyn}}$ is the above-
defined dynamical time of the disc, and $\Delta T_{\rm{impacts}}$ is the time
interval between two impacts of Sgr. As shown in Figure 1 (bottom), the ratio
$T_{\rm{dyn}}/\Delta T_{\rm{impacts}}$ demarcates four different regimes:

  * •

Regime 0 or pre-interaction regime: $T_{\rm{dyn}}/\Delta T_{\rm{impacts}}\sim
0$, as $\Delta T_{\rm{impacts}}\sim\infty$. The disc effectively evolves in
isolation.

  * •

Regime 1: $0<T_{\rm{dyn}}/\Delta T_{\rm{impacts}}<1$. In this Regime,
interaction time-scales are long enough to give the disc the possibility to
relax between one and the other plane-crossing.

  * •

Regime 2: $T_{\rm{dyn}}/\Delta T_{\rm{impacts}}\sim 1$. Since the disc
dynamical time is of the same order of the time interval between two impacts,
the disc doesn’t have time to fully relax after an interaction.

  * •

Regime 3: $T_{\rm{dyn}}/\Delta T_{\rm{impacts}}>1$. This is the most perturbed
regime, where the repeated impacts on short time-scales in the inner disc
cause a strong and complicated response in the Galactic disc.

As shown in Figure 1 (bottom), the interaction regime strongly depends on
radius R. As a result, for a given time, different regimes may coexist at
different radii.

While this simulation is not an exact reconstruction of the MW-Sgr
interaction, it is useful to consider what point in the simulation is most
similar to the present-day Milky Way-Sgr system. Based on current data (Ibata
et al., 1994; Johnston et al., 2005; Law & Majewski, 2010; Peñarrubia et al.,
2010), Sgr’s position in the Milky Way is just below the Galactic plane, at
about $\sim 16$ kpc from the Galactic center. Considering the L2 model, the
most similar geometry is at the 4th and 5th plane-crossings (i.e. R$\sim 10$
kpc), which occurs between Regime 2 and 3.

Figure 3: Face-on view of the mean vertical distortion of the Galactic disc
during the first interaction with Sgr (first column), accompanied by the
corresponding Fourier decomposition into the $m=0$ (second column), $m=1$
(third column) and $m=2$ (fourth column) terms. The Galaxy rotates clockwise.
The snapshot at t=2.35 Gyr has been marked by a black star, to highlight that
an interaction with Sgr is taking place (the pericenter and plane-crossing
occur, respectively, at 2.33 and 2.36 Gyr). The Fourier decomposition shows
that: (i) Sgr excites ring-like vertical distortions (second column), which
propagate from the inner to the outer parts of the disc; (ii) The $m=1$
distortion moves in a prograde direction with respect to Galactic rotation;
(iii) The $m=2$ term is excited by the interaction with Sgr. An animated
version of this figure is available for the entire simulation.  Figure 4: For
a given time t and Galactocentric radius R, the Fourier term with the largest
amplitude based on our analysis. Dashed vertical lines indicate plane-
crossings of Sgr. The position of Sgr during a plane-crossing is indicated by
a white diamond if $R_{\rm{Sgr}}<25$ kpc. Otherwise, a black arrow is drawn at
the top of the figure.

Figure 5: Amplitude of the Fourier vertical terms m=0 (upper panel), m=1
(middle panel) and m=2 (lower panel). Dashed vertical lines, arrows and white
diamonds are as in Figure 4. Figure 6: Amplitude of the vertical Fourier m=0
term. Dashed vertical lines, arrows and white diamonds are as in Figure 4.

##  4 Vertical response of the disc

###  4.1 Qualitative overview

Here, we qualitatively discuss the response of the disc to the repeated
impacts of Sgr. Although it is not the main objective of this paper, it is
instructive to observe the evolution of the disc stellar density (Figure 2,
first row). The disc is initially smooth (first panel), and develops a grand
design spiral structure after the first interaction (second panel). After
three disc crossings, spiral arms are not as well-defined as before, and
present small-scale discontinuities in the outer parts (third panel). Finally,
after eight disc crossings, the spiral arms are disrupted and a central bar is
formed (fourth panel).

Different interaction regimes clearly manifest themselves also in the vertical
dimension (Figure 2, second row). The disc initially exhibit a mild S-shape
warp (first panel), just before the first disc crossing. After the first
interaction with Sgr, the distortion of the disc remains relatively smooth,
but spiral waves are excited (second panel, see following Sections for a
deeper analysis). Such waves can be clearly seen for R < 15 kpc, are trailing,
but do not spatially coincide with the spiral arms observed in the density. As
plane crossings become more frequent and Sgr gets closer to the disc, the
vertical distortion becomes more prominent (third panel), until vertical
features loose their coherence over large portions of the Galactic disc
(fourth panel).

A similar behaviour is observed in the mean vertical motions $V_{Z}$ (Figure
2, third row), which appear to be related to the vertical distortion in Z (see
Figure 2, second row, and later discussion on the connection between Z and
$V_{Z}$).

A qualitative impression from Figure 2 already gives a taste of the response
of the disc to different interaction regimes, but a quantitative analysis is
needed to understand the physical mechanisms governing the evolution of the
disc under the repeated interactions with a satellite. In the following
Sections, we will concentrate on the vertical dimension, quantitatively
analyze its temporal evolution and discuss its connection with Sgr’s orbit.

###  4.2 Fourier analysis of vertical oscillations

In this Section, we quantify the vertical response of the disc by means of a
Fourier analysis. For a given ring at Galactocentric radius R, we express the
mean Z as a sum of infinite Fourier terms terms in Galactic azimuth $\phi$:

$Z(R,\phi)=A_{0}(R)+\sum_{m=1}^{\infty}A_{m}(R)\cos(m\phi-\Phi_{m}(R))\quad,$
(1)

where $A_{0}(R)$ is the amplitude of the Fourier term $m=0$ at a given R,
which is a constant with respect to $\phi$, while $A_{m}(R)$ and $\Phi_{m}(R)$
are, respectively, the amplitude and phase of the vertical terms $m\geq 1$.
The azimuthal angle $\phi$ is taken as positive in the direction of Galactic
rotation. Rings in R are created in such a way that they contain the same
number of particles (40 000), implying that they are more concentrated in the
inner regions, and less dense in the disc outskirts. This causes the most
external ring to be at $\simeq 23$ kpc, beyond which the requested number of
particles cannot be reached.

In this paper, we follow the method used by Chequers et al. (2018); Chequers &
Widrow (2017) to study m=1 bending waves in an isolated Milky Way-like disc
and under the continuous bombardment of sub-halos. However, our approach is
somewhat complementary to their analysis, as we analyse the response of the
disc to a single massive Srg-like perturber, with the specific goal of
identifying typical behaviours in different dynamical regimes of interaction.

Figure 3 illustrates the power of Fourier decomposition. For a given radius R,
every term $m$ can be summarized by either one (for $m=0$) or two (for $m\geq
1$) numbers, indicating, respectively, the amplitude, or amplitude and phase.
This makes the approach particularly convenient, because we can reduce to some
extent the degree of complexity observed in the disc vertical response (first
column) by applying a Fourier decomposition, and analyzing the temporal
evolution of the amplitudes and phases of the single terms separately (second,
third and fourth column). As will become apparent in the following, this
approach allows us to identify several features, which otherwise would have
been buried within the complexity of the disc response.

For practical reasons, we first aim to understand whether some Fourier terms
are more relevant than others, so that we can devote our attention to the
dominant ones, and observe their evolution over time. Figure 4 shows which
Fourier vertical term is dominant at a given time and a given distance from
the Galactic center. Vertical dashed lines show Sgr’s plane crossings. The
plot summarises in a compact manner the evolution of the disc: a vertical
slice corresponds to a single snapshot of the simulation, while a horizontal
slice shows the temporal evolution at a given annulus at R in the Galactic
disc. We can clearly see that the most relevant terms are $m=0,1$ and 2. In
contrast, the mode $m=3$ makes a minor contribution, and higher terms are
completely absent. Based on the above considerations, in the rest of the paper
we will focus on the first three terms, $m=0,1$ and 2.

####  4.2.1 Overview of the responses in different regimes

The separate evolution of the $m=0,1$ and 2 amplitudes is shown in Figure 5.
As we can see, in general, interactions with Sgr are typically associated with
a significant increase in the amplitude, for all three terms. Then, the
amplitude oscillates a few times, before going back to a semi-quiescent phase.
The oscillatory behaviour induced by the interaction with Sgr is more
significant and lasts longer in the outer parts of the disc than in the inner
regions.

In the last part of the simulation, the dynamical time-scales of the
satellite’s orbit are not long enough for the outer disc to recover from each
interaction. This effects manifests itself as a rapid increase after one
interaction, followed by the above-mentioned oscillatory behaviour, but with a
typically very short (Regime 2) or event absent (Regime 3) quiescent phase
before another impact occurs.

It is also worth noting that not all plane-crossings have the same impact on
the disc’s behaviour. For example, during the plane-crossings at $\approx 4.9$
and $6.1$ Gyr, Sgr is far from pericenter (Figure 1), and the disc reaction is
less relevant (or even irrelevant) compared to the other cases.

Finally, we note that the repeated impacts of Sgr on the Galactic disc cause a
general amplitude increase of the Fourier vertical terms m=0,1 and 2
throughout the simulation, implying typically stronger vertical distortions in
the more perturbed regimes (Regime 2 and 3) than in the quiescent phases.

####  4.2.2 Fourier vertical term m=0

Unlike the other Fourier terms (where $A_{m}(R)>0$ by definition), the
amplitude of $m=0$ can be either positive or negative. Figure 6 repeats the
top panel of Figure 5, but this time color-coded by sign. In the outer parts
of the disc, oscillations are clearly visible. Late regimes (Regime 2 and 3)
exhibit amplitudes up to 2 order of magnitudes larger than the early ones
(Regime 0 and 1).

In response to an interaction, the $m=0$ term manifests itself as a series of
ring-like vertical perturbation of opposite sign, typically propagating from
the inner to the outer disc, similar to a stone causing ripples to propagate
out in all directions across the surface of a pond. Ring-like vertical
distortions are propagating faster close to an interaction, and become much
slower or almost static in quiescent phases. If the sign of $A_{0}$ is
constant for all radii and the amplitude increases with R (e.g. vertical
monochromatic features at $\approx$ 4.6 Gyr or 5.7 Gyr), the $m=0$ mode forms
an U-shaped bending mode (e.g. Reshetnikov & Combes, 1998). This can happen in
quiescent phases or in the immediate proximity of an interaction (before
evolving into the above-mentioned outgoing vertical perturbations).

A visual representation of the $m=0$ temporal evolution can be found in Figure
3, 17 and 18.

Figure 7: Evolution of the phase angle for the m=1 Fourier term. Dashed
vertical lines, arrows and white diamonds are as in Figure 4. A zoom of the
last region is given in the Appendix (Figure 19).

####  4.2.3 Fourier vertical term m=1

The m=1 term corresponds to a classical S-shape warp (e.g. Chen et al., 2019;
Skowron et al., 2019) and might be the dominant Fourier term in the Milky Way
(Levine et al., 2006). Following the approach of Chequers et al. (2018),
Figure 7 shows the the temporal evolution of the m=1 phase angle. Such a plot
can be used to determine the morphology of m=1 bending waves in the disc.
Provided that the Galactic azimuth $\phi>0$ is positive in the direction of
Galactic rotation, one can dissect the plot horizontally, and observe the
variation of the phase angle for a fixed radius R as a function of time.
Hence, one can deduce if the phase angle (i.e. line-of-nodes of the warp) is
moving in a prograde (towards $\phi>0$) or retrograde (towards $\phi<0$)
direction with respect to Galactic rotation. Alternatively, one can dissect
the plot vertically, and determine whether the line-of-nodes of the warp is
straight (i.e. constant with R), leading (i.e. $\phi$ increases with
increasing radius R) or trailing (i.e. $\phi$ decreases with increasing radius
R).

After an initial relaxation time due to the non-equilibrium initial
conditions, a clear coherent pattern manifests itself between $t\approx 0.9$
and $\approx 2.1$ Gyr. At $R\approx 15$ kpc, the phase angle progressively
changes from $\approx-45$ deg (dark blue color) to $-170$ deg (light blue). A
similar behaviour is observed at all Galactic radii, but especially in the
outer parts, with the phase angle slowly moving towards $\phi<0$, thus
opposite to Galactic rotation. Moreover, for a given time, we can take a
vertical slice of the plot, and observe that the phase angle varies towards
$\phi>0$ for increasing Galactic radii R, forming a leading spiral. For
example, at $t\approx 2$ Gyr the phase angle is $\approx-180,-140$ and $-100$
deg, respectively, at $\approx 12,16$ and 20 kpc. It is worth noting that this
pattern is characterized by a long time-scale, which causes the line-of-nodes
to vary slowly, with a period much larger than 1 Gyr. This behaviour has been
extensively discussed in the literature (Shen & Sellwood, 2006; Gómez et al.,
2016; Chequers et al., 2018).

The first plane-crossing of Sgr brings the above-described pattern to an end.
The interaction with Sgr resets the phase angle at all radii, and causes the
outer parts of the disc ($R\gtrapprox 8$ kpc) to precess (term used in
literature to indicate the variation of the line of nodes) in a fast and
almost rigid way in the same (prograde) direction with respect to Galactic
rotation. At $15$ kpc, the prograde behaviour is maintained for about $200$
Myr. Right after, abrupt variations of the phase angle take place,
corresponding to the oscillatory behaviour observed in Section 4.2.1 (i.e.
disc flapping). Visually, $m=1$ prograde distortions with a different
orientation (i.e. phase angle) with respect to those in the outer parts appear
to originate from the inner parts and move to the outer parts of the disc. In
contrast, the inner parts of the disc ($R\lessapprox$ 8 kpc) exhibit a
retrograde precession, which is progressively restored also in the outer parts
of the disc, as soon as the perturbation caused by the interaction is damped
over the course of time. Therefore, Regime 1 is characterized by a slow
recovery from the first interaction with Sgr, where the disc progressively
goes back to a slow, retrograde, differential pattern speed. Also in this
case, the line-of-nodes progressively winds up into a leading spiral.

The second plane-crossing of Sgr (which, in this case, approximately coincides
with its second pericentric passage) marks the end of Regime 1. From this
point on, time scales become faster, and features are often not coherent over
all Galactic radii. It is therefore convenient to zoom in the last region, as
shown in Figure 20 (top).

At about 4.5 Gyr, after Sgr’ second plane-crossing, most of the disc ($R>5$
kpc) exhibit a fast, prograde and almost rigid precession, in agreement with
what was observed during the first interaction. Between 10 and 20 kpc, the
average period (i.e. time required for the line-of-nodes to make one complete
rotation around the Galactic center) is of about 300 to 400 Myr, corresponding
to an angular speed of about $20$ km/s/kpc. (This information can be
approximately obtained from Figure 20 by estimating the time difference
between two recurrences of the same phase angle. However, pattern speeds will
be more extensively discussed in the following Sections.)

Between 5.5 and 6.1 Gyr, the outer parts of the disc respond to the impact of
Sgr in a similar way to the interactions described before, while the inner
parts seem to develop the slow, retrograde pattern observed in Regimes 0 and 1
after about 5.7 Gyr. This situation represents an illustrative example of how
different regimes may co-exists at the same time, due to different dynamical
time-scales in the Galactic disc. Indeed, here the time interval between two
impacts is of the same order of the disc dynamical times in the outer disc
(Regime 2), but still larger than those in the inner disc (Regime 1).

After 6.1 Gyr, plane-crossings occur so frequently that the disc is
continuosly subjected to external perturbations. As a consequence, complicated
and discontinous patterns arise in the disc, being clearly visible also in the
inner regions, in which Sgr is progressively being absorbed. In the outer disc
($\approx$ 10 - 20 kpc), a coherent structure survives from about 6.2 to 6.7
Gyr, showing a relatively fast, prograde and nearly rigid precession,
completing half of the disc ($180$ deg) in about 300 Myr (roughly
corresponding to an angular speed of $\approx$ 10 km/s/kpc). Similarly, other
features may survive for a few hundreds of Myr, but they do not remain
coherent for more than a few kpc in Galactic radius. Therefore, the recurrent
impacts on a short time-scale, together with the proximity of the satellite,
make Regime 3 the most perturbed regime analyzed in this work, dominated by
complex disc patterns, which do not remain coherent for very long, both in
space and time.

A visual representation of the $m=1$ slow retrograde pattern, typical of
quiescent phases, can be seen in Figure 17. The $m=1$ prograde pattern, which
occurs during an interaction, can be visualised in Figure 3 and 18.

Figure 8: Evolution of the phase angle for the m=2 Fourier term. Dashed
vertical lines, arrows and white diamonds are as in Figure 4. A zoom of the
last region is given in the Appendix (Figure 19).

Figure 9: Difference between the phase angle of the m=1 and m=2 Fourier terms
at different Galactocentric radii R. Dashed vertical lines as in Figure 4.
Close to interactions, the two terms are moving at the same speed (i.e. where
the difference is constant), illustrating the coupling of the Fourier terms.
In the bottom panel, the phase angle inside 12 kpc are represented as open
symbols, to avoid crowding.

####  4.2.4 Fourier vertical term m=2

The m=2 term corresponds to a bi-symmetric vertical distortion of the Galactic
disc. In general, its temporal evolution presents some relevant differences
from the previously discussed $m=1$. For example, $m=2$ is completely absent
in the pre-interaction regime, and gets excited by the first plane crossing of
Sgr (at $\approx$ 2.3 Gyr), as shown by Figure 5.

The temporal evolution of the $m=2$ phase angle is shown in Figure 8 (see also
Figure 20, which is zoomed in the most perturbed region). As we can see, the
$m=2$ term exhibits a very regular behaviour. Even in the latest regimes,
where Sgr’s impacts become more frequent and the disc is strongly perturbed,
the $m=2$ pattern can be roughly recognised. For comparison, the m=1 term
appears to be much more fragile and easily disrupted.

As will be better explained in the spectral analysis, the $m=2$ bending waves
typically move in a prograde direction (also in quiescent phases, unlike
$m=1$). Their angular speeds show a very mild variation with radius, which
causes the bending waves to wrap up very slowly. Starting from an initially
long radial wavelength ($>25$ kpc), the phase of the m=2 term progressively
winds up, reaching a radial wavelength of $\approx$ 8 kpc almost 2 Gyr after
the latest interaction. Therefore, the pitch angle of the m=2 vertical wave
can, in principle, be used as a clock to date the most recent interaction.

During a plane-crossing, if that also coincides with a pericentric passage,
the phase-angle of the m=2 term gets reset at all radii, and, intially
straight, moves in a prograde direction, similarly to what observed for the
$m=1$ term. Figure 9 shows the difference between the phase angle of the $m=1$
and the $m=2$ term ($\phi_{1}$ and $\phi_{2}$, respectively) at different
times. During the first two pericenters, the $m=1$ and the $m=2$ term seem to
be somewhat aligned, moving at approximately the same speed (i.e. their
difference is roughly constant). Such behaviour might be due to a coupling
between the Fourier terms, which occurs as a response to the interaction with
Sgr. However, this effect becomes less clear in the latest and more perturbed
regimes, where the presence of Sgr inside the disc does not allow coherent
features to survive for more than a few kpc in Galactic radius. In contrast,
very regular patterns can be observed in quiescent regions (e.g. $3<$ t $<4.4$
Gyr), where the two Fourier terms revolve at their own typical speed, with the
$m=1$ being slightly retrograde (see above) and the $m=2$ typically prograde.
This difference is causing the difference $\phi_{1}-\phi_{2}$ to be decreasing
with time.

Figure 10: Vertical spectrograms of the Fourier m=1 vertical term. Each panel
covers a temporal bin of 1 Gyr in the simulation. Panels marked by a black
star are those where one (or more) interaction(s) with Sgr took place. The
solid curve shows corotation. The dashed curves show the vertical resonances
at $\Omega(R)\pm\nu_{z}(R)$, where $\Omega(R)$ is Galactic rotation and
$\nu_{z}(R)$ is the vertical frequency (see Appendix). The dashed-dotted
curves show the vertical resonances at $\Omega(R)\pm 2\nu_{z}(R)$,  Figure 11:
Same as Figure 10, but for the Fourier m=2 vertical term.

####  4.2.5 Spectral analysis

Originally developed for density waves (Sellwood & Athanassoula, 1986; Masset
& Tagger, 1997; Quillen et al., 2011; D’Onghia & L. Aguerri, 2020), spectral
analysis has been recently extended to bending waves (Chequers et al., 2018;
Chequers & Widrow, 2017), allowing us to describe them in terms of frequency
power spectrum as a function of Galactocentric radius. A full description of
vertical spectrograms can be found in Chequers et al. (2018); Chequers &
Widrow (2017), and technical details are given in Appendix C. In this paper,
we will use vertical spectrograms to identify the most relevant features over
a given time baseline. Specifically, we will be interested in determining the
angular speed (i.e. temporal variation of the phase angle $\Phi_{m}(R)$, see
Equation 1) of the most relevant Fourier terms $m\geq 1$ in different regimes
of interaction with Sgr.

Figure 10 consists of six vertical spectrograms, showing the frequency
spectrum for $m=1$ bending waves in different periods of the simulation. To
relate the vertical response of the disc to the satellite’s orbit, we chose a
temporal baseline of 1 Gyr for each spectrogram, which implies a resolution on
the vertical axis of 6.28 km/s/kpc (the temporal resolution of the simulation
is 10 Myrs). (For a sharper resolution on the vertical axis, implying a longer
baseline, see Appendix C.) One of the most obvious features of Figure 10 is
that the power reaches typically larger values (i.e. stronger bending waves)
in the spectrograms where one (or more) interaction(s) with Sgr occurred
(panels marked by a black star), whereas it remains relatively low in
quiescent regions (panels A and C).

It is also instructive to look at where the power is concentrated at different
stages of the Milky Way-satellite interaction. In quiescent regions, most of
the power is concentrated close to angular speed 0 (panel A), or slightly
retrograde values (panel C), in the outer regions of the disc. This
approximately aligns with the vertical resonance $\Omega-\nu_{z}$, where
$\Omega$ is the angular rotation rate of the Galaxy and $\nu_{z}$ is the
vertical epyciclic frequency. This is the so-called slow wave, and is expected
to be retrograde because, in a flattened galaxy, the vertical frequency
$\nu_{z}>\Omega$. The slow and retrograde bending wave has been already
discussed in the literature (Binney & Tremaine, 2008; Chequers & Widrow, 2017;
Chequers et al., 2018).

However, in panel B, power is located in two main regions: one close to
$\simeq$20-25 km/s/kpc, indicating the prograde precession associated with
Sgr’s impact, and the slow retrograde precession, which is visible
before/after the effects of the perturbation have relaxed. The prograde
feature disappears in panel C, where no interaction take place. Therefore,
there is a clear connection between the vertical pattern speeds in the
Galactic disc and the satellite’s orbit. This is further confirmed by panels
D, E and F, where multiple interactions are present, and transient prograde
features become prominent. We also observe that prograde angular speeds in
panel E and F are typically slower than those in panel B and D.

Figure 11 repeats the vertical spectrograms of Figure 10, but now for the
$m=2$ term. Considerations similar to the ones above can be done here, with
the difference that: (i) the m=2 term is completely absent before the first
interaction (panel A), and (ii) the resonance curves are different, as the
angular speed $\omega=m\Omega\pm\nu_{z}$, with $\omega=m\Omega_{p}$, where
$\Omega_{p}$ is the pattern speed of the bending wave. In the previous
sections, we found that $m=2$ bending waves present a very regular behaviour,
characterized by a prograde angular speed. Indeed, the spectrograms indicate
that, in the outer regions, the power aligns with the vertical resonance at
$\omega=m\Omega-\nu_{z}$ (see Bland-Hawthorn & Tepper-Garcia, 2020), which
demarcates a prograde angular speed.

Figure 12: Comparison between the phase angle of the m=2 mode for Z (red), Vz
(blue) and density (black) at t=3.3 Gyr, 5.5 Gyr and 6.5 Gyr (from left to
right).

Figure 13: Same as Figure 9, but showing the difference between the m=2
density and m=2 vertical Fourier terms.

####  4.2.6 Relating response in different dimensions

Figure 12 shows the phase angle of the $m=2$ vertical term (Z) as a function
of Galactocentric radius R, compared to other quantities that might be of
interest. The $m=2$ term in vertical velocity ($V_{z}$) appear to be typically
shifted by $90\deg$ with respect to the $m=2$ phase in Z, as expected from a
wave-like motion.

When comparing the $m=2$ vertical term (Z) and the $m=2$ spiral arms in
density, we find that, in general, there is no obvious correlation. This can
be also seen in Figure 2, where the spiral arms in density appear more tightly
wound than the ones in the vertical spiral wave. Indeed, such behaviour is
somewhat expected, as the $m=2$ density and bending waves have, respectively,
a pattern speed of $\Omega(R)-\kappa/2$ and $\Omega(R)-\nu_{z}/2$. This is in
good agreement with the recent analysis presented by Bland-Hawthorn & Tepper-
Garcia (2020).

Figure 13 shows the angular difference between the $m=2$ phase angle in
density and in Z as a function of time for different Galactocentric radii. As
we can see, the $m=2$ density and bending waves usually move at their own
typical rate, with the $m=2$ in density wave being typically faster than $m=2$
in $Z$. This is apparent when the difference $\phi_{2}$ (Density) - $\phi_{2}$
(Z) increases as a function of time. However, in proximity of an interaction,
the behaviour is different. The $m=2$ density and bending waves tend to align,
and move at a similar speed for a very short time. This behaviour might
highlight a resonance, caused by the interaction with Sgr.

Figure 14: Comparison between the angular speed obtained from a data-like
approach (red points, see text) and the actual behaviour of the disc (grey
points and shaded area) for the m=1 (left) and m=2 (right) Fourier vertical
terms. Different rows refer to different positions in the Galactic disc. The
grey points are obtained by simply subtracting the phase angles at two
subsequent snapshots (for a more sophisticated approach, see Section 4).

Figure 15: Same as Figure 14, but comparing the results of the data-like
approach using the entire disc (red points) and a limited volume of the
Galactic disc (blue points). For a given snapshot, the blue points show the
median angular speed obtained from eight different volumes (see text), while
the error bar shows the median absolute deviation (MAD) of the obtained
values.

##  5 Observing the disc response

The previous Section showed that the evolution of the Galactic disc and the
orbit of a perturbing satellite are tightly correlated, thereby leading to
various expectations for vertical pattern speeds in the Galactic disc and
other quantifiable properties, depending on the status of interaction. In
principle, the response of the disc can be inferred from observations, and
used to learn more on the interaction with the perturbing satellite(s) and the
perturbed halo, with the caveat that other internal mechanisms should be
considered as well, such as the bar buckling (with the further caveat that we
do not have an indication of when in the lifetime of the Galaxy this may have
occurred).

The evolution of the disc in N-body simulations can be deduced by comparing
subsequent snapshots at different times (as was done in the previous Section).
However, in real data, observations provide us with positions and velocities
of stars in the Galaxy only at the present time. Therefore, a method is
required to infer the actual behaviour of the disc using only the
instantaneous measurements of the 6-dimensional phase-space of stars.

In this Section, we aim to mimic the real data approach, in order to
understand whether it is possible to recover the actual instantaneous
behaviour of the Galactic disc by simply analysing positions and velocities of
stars at a given time. Certainly, the equations governing the vertical
response of the Galaxy to the recurrent impacts of a satellite galaxy should
be based on a deep dynamical investigation, which is beyond the scopes of the
present paper. Instead, here we use a very simple kinematical approach, and
test its validity and limitations. The approach currently includes the three
Fourier terms $m=0,1$ and 2 (Appendix A shows that this is a reasonable
approximation), but it can be easily extended to other Fourier terms, if
needed. We model the disc using an exponential vertical profile, whose
vertical displacement is determined as follows. For a given R, the amplitude
$h_{0}$ of the $m=0$ term is free to change with time. The two terms $m=1$ and
$m=2$ have a time-varying amplitude ($h_{1}$ and $h_{2}$, respectively) and
phase ($\phi_{1}$ and $\phi_{2}$, respectively). The above time-dependent
vertical displacement has implications on the vertical velocities of the
stars, whose mean trend can be calculated using the Collisionless Boltzmann
Equation (see Equation 6 and its derivation for the $m=1$ term in Poggio et
al., 2020):

$\begin{split}V_{z}(R,\phi)=\frac{\partial h_{0}}{\partial
t}+\Bigl{(}\frac{V_{\phi}}{R}-\frac{\partial\phi_{1}}{\partial
t}\Bigr{)}h_{1}\cos(\phi-\phi_{1})+\frac{\partial h_{1}}{\partial
t}\sin(\phi-\phi_{1})\\\
+\Bigl{(}2\frac{V_{\phi}}{R}-\frac{\partial\phi_{2}}{\partial
t}\Bigr{)}h_{2}\cos(2\phi-\phi_{2})+\frac{\partial h_{2}}{\partial
t}\sin(2\phi-\phi_{2})\quad,\end{split}$ (2)

where the kinematic parameters $\frac{\partial\phi_{1}}{\partial t}$ and
$\frac{\partial\phi_{2}}{\partial t}$ are the angular frequencies of the $m=1$
and $m=2$ terms, respectively, and $\frac{\partial h_{0}}{\partial t}$,
$\frac{\partial h_{1}}{\partial t}$, $\frac{\partial h_{2}}{\partial t}$ are
the amplitude variations of the corresponding Fourier terms. These five
parameters are, in general, non-zero, and can be statistically inferred using
positions and velocities ($R,\phi,V_{\phi},V_{z}$) of the stars in the
Galactic disc at a given time using Equation 2. To do this, in this experiment
we use a simple non-linear least squares method, assuming that the (spatial)
amplitudes and phases of the $m=0,1$ and 2 terms are known. Hereafter, we will
refer to this method as the _Instantaneous measurement_ , since the kinematic
parameters are recovered from the 6-dimensional phase-space information for
stars at a fixed time, as in the case of real data.

As anticipated above, an alternative way to determine the instantaneous
evolution of the Galactic disc is to compare the amplitudes and phases at
different snapshots. For example, one can calculate the angular speed
$\frac{\partial\phi_{1}}{\partial t}$ at given radius R by naively taking the
phase difference at two subsequent snapshots. Alternatively, one might use the
vertical spectrograms (see Section 4.2.5), which, however, provide us with the
power distribution on a relatively large time-scale (1 Gyr), due to the
simulation temporal resolution. Therefore, in this simple experiment, we use
the first approach, and name it _Differential Fourier analysis_.

Figure 14 shows the results obtained when recovering the angular speeds
$\frac{\partial\phi_{1}}{\partial t}$ and $\frac{\partial\phi_{2}}{\partial
t}$ of the $m=1$ and $m=2$ terms according to the Instantaneous measurement
and the Differential Fourier analysis at different Galactic radii. As we can
see, notwithstanding its simplicity, our data-like approach seems to produce
results in rough agreement with the Differential Fourier analysis, even in the
latest regimes, where the disc is more perturbed.

However, a realistic magnitude limited survey does not cover the entire
Galactic disc. Therefore, we perform an additional experiment, to test the
applicability of the data-like approach on a limited portion of the Galactic
disc. We place the simulated Sun at 12 kpc from the Galactic center, to
approximately reproduce a realistic warp amplitude based on the observed mean
vertical displacement at a given heliocentric distance (e.g. Momany et al.,
2006). While the Sun’s distance from the Galactic center is maintained fixed,
we place the Sun at 8 different azimuthal angles (spaced by 45o) for a given
snapshot. Then, we select only stars within 10 kpc from the Sun, to
approximately mimic a Gaia DR3-like sample of stars. For each snapshot, we end
up with 8 mock catalogs of stars, containing $\simeq 600\,000$ particles
(which implies a lower statistic than the one that can presumably be reached
with Gaia DR3). Due to the low statistics, in this second experiment we fit
for a single value of angular speed $\frac{\partial\phi_{1}}{\partial t}$ and
$\frac{\partial\phi_{2}}{\partial t}$, instead of evaluating each radial ring
separately. (This approximation is equivalent to averaging over the portion of
the disc covered by the mock dataset.) Finally, we take the median of the
results from the 8 mock catalogs, and calculate the dispersion around the
median value, and repeat the procedure for all the snapshots of interest. The
results are presented in Figure 15. As we can see, the dispersion can be
significant in some places, as expected, due to the complexity of the physical
mechanisms at work in our Galaxy, which are presumably not captured by our
simple approach. However, the median value is in good agreement with the
values obtained using the entire disc, meaning that the approach is unbiased.

##  6 Discussion and conclusion

###  6.1 Key results

Using N-body simulations from L18, we have described the interaction between
the Milky Way and Sgr in terms of four different interaction regimes (called
Regime 0,1,2,3), based on the orbital history of Sgr and the dynamical time
scales in the Galactic disc. Each regime dictates the properties, extent and
evolution of the disc perturbations. We have quantitatively characterised the
evolution of the vertical distortions of the disc under the repeated impacts
of Sgr using Fourier decomposition and spectral analysis. Some of the observed
features have been already discussed in literature, while others are new.
Finally, we have tested the validity of pattern speed measurements in a
realistic model. The main results of the paper are summarized below:

  1. (1)

The $m=0$ term manifests itself as a ring-like vertical distortion with
alternate sign propagating from the inner to the outer parts of the disc. The
radial propagation speed is faster after an interaction, and becomes much
slower (almost static) in quiescent phases.

  2. (2)

The $m=1$ term is characterized by a slow, differential, retrograde precession
in the quiescent regimes, while it moves in a relatively fast and prograde
direction in the last regimes or close to an interaction.

  3. (3)

The $m=2$ term is excited by the first interaction with Sgr, and presents a
very regular behaviour, apparent as a prograde bi-symmetric trailing spiral,
which progressively winds up with time after an interaction.

  4. (4)

In quiescent phases, the Fourier vertical terms $m=1$ and $m=2$ revolve at
their own typical rate, which is determined by vertical Lindblad resonances.
However, their angular frequency is forced to change during the interactions
with Sgr, and they both move in a prograde direction for a few hundred Myrs,
at approximately the same speed. Such coherent motion might perhaps indicate a
resonance, which is arising as a response of the disc to the interaction with
the satellite.

  5. (5)

The instantaneous evolutionary state of the Galactic disc can be parameterised
using the amplitude and phase angle variation (angular speed) of the Fourier
terms. Such parameters can be recovered in an unbiased way from stellar
positions and velocities in the Galactic disc at a given time (i.e. similar to
the real data case).

###  6.2 Implications

In this work, we have shown that the $m=0,1$ and 2 Fourier vertical terms
dominate and behave in a different way in distinct phases of the interaction
process with a satellite galaxy. Although some caveats (discussed below)
should be kept in mind, our results imply that _the kinematic parameters
describing the instantaneous vertical evolution of the Galactic disc have a
constraining power on the dynamical state of the interaction_ , when a Milky
Way-Sgr like accretion is considered. Taking advantage of current and future
datasets, which allow us to characterize the Galactic disc on an
unprecedentedly large volume, one can (i) place observational constraints on
the evolution of the perturbed Galactic disc, and (ii) use high-resolution
N-body simulations to bridge the gap between observation and theory, to
suggest a possible dynamical interpretation of the observations in the context
of a Milky Way-Sgr interaction.

The scenario analysed in this work provides a possible framework for the
interpretation of the measured warp ($m=1$) pattern speed ($\simeq 11\pm 3$
km/s/kpc, 13.6 km/s/kpc, respectively Poggio et al., 2020; Cheng et al.,
2020), which would be more consistent with the Milky Way currently being in
the second or third dynamical regime (see Figure 10E, 10F, and 14) or, in any
case, close to an interaction, rather than in a semi-quiescent phase (e.g.
Figure 10A and 10c, where no interactions have occurred in a time-interval of
$\simeq 1.5$ Gyr).

Additional measurements can potentially give further insights into the origin
of vertical disturbances and relate to satellite contribution. For example,
the pattern speed of the $m=2$ vertical term might be measured using future
datasets. Such result would be not only interesting by itself, but also
because the difference between the angular speed of the m=1 and m=2 terms
represents a complementary way of determining if an interaction has occurred
in the last few hundreds of Myrs. Additionally, given that the $m=2$ mode is
persistent and winds-up very slowly (see Fig. 8), its pitch angle, even
measured locally, may provide a good clock for measuring the time since the
last perturbation.

###  6.3 Limitations and future prospects

This work aims to make further progress towards the understanding of the
vertical response of the disc to the recent accretion history of the Galaxy.
As a starting point, we analysed the response to the repeated impacts of a
Sgr-like satellite. However, a number of caveats should be kept in mind, when
interpreting our results in their wider context. For instance, modelling Sgr
orbital history from the virial radius crossing to the present-day is
notoriously difficult, which is further exacerbated by the large uncertainties
on its parameters (e.g. progenitor mass profile, last pericentric passage,
etc.), and those of the Milky Way’s own potential. Hence, instead of aiming to
reproduce the _exact_ orbital history of Sgr and the _exact_ vertical response
of the disc, we have rather concentrated ourselves on the exploration of the
physical spatial reaction to (and the relative time-scales of) an interaction
between a Milky Way-like galaxy and a Sgr-like satellite. In addition to this,
it should be noted that the Magellanic Clouds also play a role in generating
vertical pattern speeds in the Galactic disc (e.g. Weinberg & Blitz, 2006;
Laporte et al., 2018a). Its effect in combination with Sgr was already
explored in Laporte et al. (2018b) showing that the recent interaction with
the LMC modulates the response due to Sgr which is particularly noticeable in
the outer disc ($\geq 16$ kpc). In addition to this, we note that the
simulation does not include a gaseous component, which might contribute to the
bending wave damping. A deeper exploration is left for future work. Finally,
it should be mentioned that our work did not consider the effects of the bar
buckling; whether it has a significant impact or not on the outer disc should
be further explored. Future works comparing observational constrains on the
disc vertical evolution with possibly more sophisticated dynamical models will
certainly reveal further information on the recent accretion history of the
Milky Way.

## Acknowledgements

EP, CL and ED thank Lawrence Widrow for useful discussions. EP acknowledges
financial support from the Center of Computational Astrophysics, Flatiron
Institute (New York, USA), which enabled a long-term visit in 2019 to develop
this project. CL thanks the Observatoire de la Côte d’Azur for hosting him on
a long-term visit. This work used the Extreme Science and Engineering
Discovery Environment (XSEDE), which is supported by National Science
Foundation grant number OCI-1053575. This work was supported in part by World
Premier International Research Center Initiative (WPI Initiative), MEXT,
Japan. EP acknowledges the partial support of ASI through contract ASI
2018-24-HH.0, Gaia Mission - The Italian Participation to DPAC. CL
acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 852839). KVJ and DF were supported by NSF grant AST-1715582.

## References

  * Antoja et al. (2018) Antoja T., et al., 2018, Nature, 561, 360
  * Aumer & White (2013) Aumer M., White S. D. M., 2013, MNRAS, 428, 1055
  * Bailin (2003) Bailin J., 2003, ApJ, 583, L79
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Bergemann et al. (2018) Bergemann M., et al., 2018, Nature, 555, 334
  * Binney (2012) Binney J., 2012, MNRAS, 426, 1324
  * Binney & Schönrich (2018) Binney J., Schönrich R., 2018, MNRAS, 481, 1501
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Bland-Hawthorn & Tepper-Garcia (2020) Bland-Hawthorn J., Tepper-Garcia T., 2020, arXiv e-prints, p. arXiv:2009.02434
  * Bland-Hawthorn et al. (2019) Bland-Hawthorn J., et al., 2019, MNRAS, 486, 1167
  * Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
  * Burke (1957) Burke B. F., 1957, AJ, 62, 90
  * Carlin et al. (2013) Carlin J. L., et al., 2013, ApJ, 777, L5
  * Carrillo et al. (2019) Carrillo I., et al., 2019, MNRAS, 490, 797
  * Chen et al. (2019) Chen X., Wang S., Deng L., de Grijs R., Liu C., Tian H., 2019, Nature Astronomy, 3, 320
  * Cheng et al. (2020) Cheng X., et al., 2020, arXiv e-prints, p. arXiv:2010.10398
  * Chequers & Widrow (2017) Chequers M. H., Widrow L. M., 2017, MNRAS, 472, 2751
  * Chequers et al. (2018) Chequers M. H., Widrow L. M., Darling K., 2018, MNRAS, 480, 4244
  * D’Onghia & L. Aguerri (2020) D’Onghia E., L. Aguerri J. A., 2020, ApJ, 890, 117
  * D’Onghia et al. (2010) D’Onghia E., Vogelsberger M., Faucher-Giguere C.-A., Hernquist L., 2010, ApJ, 725, 353
  * D’Onghia et al. (2013) D’Onghia E., Vogelsberger M., Hernquist L., 2013, ApJ, 766, 34
  * D’Onghia et al. (2016) D’Onghia E., Madau P., Vera-Ciro C., Quillen A., Hernquist L., 2016, ApJ, 823, 4
  * Deason et al. (2018) Deason A. J., Belokurov V., Koposov S. E., 2018, MNRAS, 473, 2428
  * Debattista & Sellwood (1999) Debattista V. P., Sellwood J. A., 1999, ApJ, 513, L107
  * Dubinski & Kuijken (1995) Dubinski J., Kuijken K., 1995, ApJ, 442, 492
  * Font et al. (2001) Font A. S., Navarro J. F., Stadel J., Quinn T., 2001, ApJ, 563, L1
  * Friske & Schönrich (2019) Friske J. K. S., Schönrich R., 2019, MNRAS, 490, 5414
  * Gaia Collaboration et al. (2018a) Gaia Collaboration et al., 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2018b) Gaia Collaboration et al., 2018b, A&A, 616, A11
  * Garavito-Camargo et al. (2019) Garavito-Camargo N., Besla G., Laporte C. F. P., Johnston K. V., Gómez F. A., Watkins L. L., 2019, ApJ, 884, 51
  * Gauthier et al. (2006) Gauthier J.-R., Dubinski J., Widrow L. M., 2006, ApJ, 653, 1180
  * Gibbons et al. (2017) Gibbons S. L. J., Belokurov V., Evans N. W., 2017, MNRAS, 464, 794
  * Gómez et al. (2013) Gómez F. A., Minchev I., O’Shea B. W., Beers T. C., Bullock J. S., Purcell C. W., 2013, MNRAS, 429, 159
  * Gómez et al. (2016) Gómez F. A., White S. D. M., Marinacci F., Slater C. T., Grand R. J. J., Springel V., Pakmor R., 2016, MNRAS, 456, 2779
  * Gómez et al. (2017) Gómez F. A., White S. D. M., Grand R. J. J., Marinacci F., Springel V., Pakmor R., 2017, MNRAS, 465, 3446
  * Haines et al. (2019) Haines T., D’Onghia E., Famaey B., Laporte C., Hernquist L., 2019, ApJ, 879, L15
  * Hayes et al. (2018) Hayes C. R., et al., 2018, ApJ, 859, L8
  * Haywood et al. (2018) Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113
  * Helmi (2020) Helmi A., 2020, arXiv e-prints, p. arXiv:2002.04340
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Ibata & Razoumov (1998) Ibata R. A., Razoumov A. O., 1998, A&A, 336, 130
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Johnston et al. (2005) Johnston K. V., Law D. R., Majewski S. R., 2005, ApJ, 619, 800
  * Kallivayalil et al. (2013) Kallivayalil N., van der Marel R. P., Besla G., Anderson J., Alcock C., 2013, ApJ, 764, 161
  * Kazantzidis et al. (2009) Kazantzidis S., Zentner A. R., Kravtsov A. V., Bullock J. S., Debattista V. P., 2009, ApJ, 700, 1896
  * Kerr (1957) Kerr F. J., 1957, AJ, 62, 93
  * Khoperskov et al. (2019) Khoperskov S., Di Matteo P., Gerhard O., Katz D., Haywood M., Combes F., Berczik P., Gomez A., 2019, A&A, 622, L6
  * Laporte et al. (2018a) Laporte C. F. P., Gómez F. A., Besla G., Johnston K. V., Garavito-Camargo N., 2018a, MNRAS, 473, 1218
  * Laporte et al. (2018b) Laporte C. F. P., Johnston K. V., Gómez F. A., Garavito-Camargo N., Besla G., 2018b, MNRAS, 481, 286
  * Laporte et al. (2019a) Laporte C. F. P., Johnston K. V., Tzanidakis A., 2019a, MNRAS, 483, 1427
  * Laporte et al. (2019b) Laporte C. F. P., Minchev I., Johnston K. V., Gómez F. A., 2019b, MNRAS, 485, 3134
  * Laporte et al. (2020) Laporte C. F. P., Belokurov V., Koposov S. E., Smith M. C., Hill V., 2020, MNRAS, 492, L61
  * Law & Majewski (2010) Law D. R., Majewski S. R., 2010, ApJ, 714, 229
  * Levine et al. (2006) Levine E. S., Blitz L., Heiles C., 2006, ApJ, 643, 881
  * Li et al. (2017) Li T. S., et al., 2017, ApJ, 844, 74
  * López-Corredoira et al. (2002) López-Corredoira M., Betancort-Rijo J., Beckman J. E., 2002, A&A, 386, 169
  * López-Corredoira et al. (2020) López-Corredoira M., Garzón F., Wang H. F., Sylos Labini F., Nagy R., Chrobáková Ž., Chang J., Villarroel B., 2020, A&A, 634, A66
  * Masset & Tagger (1997) Masset F., Tagger M., 1997, A&A, 322, 442
  * McMillan (2011) McMillan P. J., 2011, MNRAS, 414, 2446
  * Meza et al. (2005) Meza A., Navarro J. F., Abadi M. G., Steinmetz M., 2005, MNRAS, 359, 93
  * Momany et al. (2006) Momany Y., Zaggia S., Gilmore G., Piotto G., Carraro G., Bedin L. R., de Angeli F., 2006, A&A, 451, 515
  * Navarro et al. (2011) Navarro J. F., Abadi M. G., Venn K. A., Freeman K. C., Anguiano B., 2011, MNRAS, 412, 1203
  * Nelson & Tremaine (1995) Nelson R. W., Tremaine S., 1995, MNRAS, 275, 897
  * Newberg et al. (2002) Newberg H. J., et al., 2002, ApJ, 569, 245
  * Peñarrubia et al. (2010) Peñarrubia J., Belokurov V., Evans N. W., Martínez-Delgado D., Gilmore G., Irwin M., Niederste-Ostholt M., Zucker D. B., 2010, MNRAS, 408, L26
  * Poggio et al. (2018) Poggio E., et al., 2018, MNRAS, 481, L21
  * Poggio et al. (2020) Poggio E., Drimmel R., Andrae R., Bailer-Jones C. A. L., Fouesneau M., Lattanzi M. G., Smart R. L., Spagna A., 2020, Nature Astronomy, 4, 590
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14
  * Price-Whelan et al. (2015) Price-Whelan A. M., Johnston K. V., Sheffield A. A., Laporte C. F. P., Sesar B., 2015, MNRAS, 452, 676
  * Purcell et al. (2011) Purcell C. W., Bullock J. S., Tollerud E. J., Rocha M., Chakrabarti S., 2011, Nature, 477, 301
  * Quillen et al. (2011) Quillen A. C., Dougherty J., Bagley M. B., Minchev I., Comparetta J., 2011, MNRAS, 417, 762
  * Ramos et al. (2018) Ramos P., Antoja T., Figueras F., 2018, A&A, 619, A72
  * Recio-Blanco (2018) Recio-Blanco A., 2018, A&A, 620, A194
  * Reshetnikov & Combes (1998) Reshetnikov V., Combes F., 1998, A&A, 337, 9
  * Rocha-Pinto et al. (2004) Rocha-Pinto H. J., Majewski S. R., Skrutskie M. F., Crane J. D., Patterson R. J., 2004, ApJ, 615, 732
  * Romero-Gómez et al. (2019) Romero-Gómez M., Mateu C., Aguilar L., Figueras F., Castro-Ginard A., 2019, A&A, 627, A150
  * Ruiz-Lara et al. (2020) Ruiz-Lara T., Gallart C., Bernard E. J., Cassisi S., 2020, Nature Astronomy,
  * Sanders & Binney (2015) Sanders J. L., Binney J., 2015, MNRAS, 447, 2479
  * Sanders & Binney (2016) Sanders J. L., Binney J., 2016, MNRAS, 457, 2107
  * Schönrich & Dehnen (2018) Schönrich R., Dehnen W., 2018, MNRAS, 478, 3809
  * Sellwood & Athanassoula (1986) Sellwood J. A., Athanassoula E., 1986, MNRAS, 221, 195
  * Semczuk et al. (2020) Semczuk M., Łokas E. L., D’Onghia E., Athanassoula E., Debattista V. P., Hernquist L., 2020, MNRAS,
  * Sheffield et al. (2018) Sheffield A. A., Price-Whelan A. M., Tzanidakis A., Johnston K. V., Laporte C. F. P., Sesar B., 2018, ApJ, 854, 47
  * Shen & Sellwood (2006) Shen J., Sellwood J. A., 2006, MNRAS, 370, 2
  * Skowron et al. (2019) Skowron D. M., et al., 2019, Science, 365, 478
  * Slater et al. (2014) Slater C. T., et al., 2014, ApJ, 791, 9
  * Sparke & Casertano (1988) Sparke L. S., Casertano S., 1988, MNRAS, 234, 873
  * Springel et al. (2008) Springel V., et al., 2008, MNRAS, 391, 1685
  * Toomre & Toomre (1972) Toomre A., Toomre J., 1972, ApJ, 178, 623
  * Vasiliev (2018) Vasiliev E., 2018, AGAMA: Action-based galaxy modeling framework (ascl:1805.008)
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 482, 1525
  * Velazquez & White (1999) Velazquez H., White S. D. M., 1999, MNRAS, 304, 254
  * Weinberg (1998) Weinberg M. D., 1998, MNRAS, 299, 499
  * Weinberg & Blitz (2006) Weinberg M. D., Blitz L., 2006, ApJ, 641, L33
  * Widrow et al. (2012) Widrow L. M., Gardner S., Yanny B., Dodelson S., Chen H.-Y., 2012, ApJ, 750, L41
  * Widrow et al. (2014) Widrow L. M., Barber J., Chequers M. H., Cheng E., 2014, MNRAS, 440, 1971
  * Williams et al. (2013) Williams M. E. K., et al., 2013, MNRAS, 436, 101
  * Xu et al. (2015) Xu Y., Newberg H. J., Carlin J. L., Liu C., Deng L., Li J., Schönrich R., Yanny B., 2015, ApJ, 801, 105
  * Younger et al. (2008) Younger J. D., Besla G., Cox T. J., Hernquist L., Robertson B., Willman B., 2008, ApJ, 676, L21
  * Yurin & Springel (2014) Yurin D., Springel V., 2014, MNRAS, 444, 62
  * de Boer et al. (2018) de Boer T. J. L., Belokurov V., Koposov S. E., 2018, MNRAS, 473, 647

##  Appendix A Visualizing the Fourier decomposition

In this Section, we give further details on the Fourier decomposition
performed in Section 4. As discussed in the text, in this paper we mainly
focus on the Fourier terms m=0, m=1, m=2 for practical reasons. Those terms
are found to be the most relevant ones (see main text), and their combination
is expected to reproduce, to a first approximation, the vertical structure of
the Galactic disc. Indeed, Figure 16 shows that the main features observed in
the simulation (first column) are already approximated in an acceptable way by
the crude sum of the first three terms (second column), although the
approximation is obviously more accurate when higher order terms are included
(third column). In this paper, higher order terms are not investigated for
practical reasons, and we leave their exploration to future works.

To complement Figure 3, which was dedicated to the first interaction with Sgr,
we here present Figure 17 and 18, showing, respectively, the evolution of the
three first Fourier terms in a quiescent phase after the first interaction
(Regime 1), and after an interaction in the most perturbed regimes (Regime 2).
As we can see, m=0 looks almost static in the quiescent phase (Figure 17),
while manifests itself as fast, outgoing, ring-like vertical distortions close
to an interaction (Figure 3 and 18). The m=1 term moves in a prograde
direction close to an interaction (Figure 3 and 18), but moves in a slow and
retrograde fashion in quiescent phases (Figure 3), creating a tightly wound
feature. After an interaction, the m=2 term progressively winds up as well,
although in Figure 18) features are more discontinuos, because of the direct
influence of Sgr (which is in the Galactic plane (z=0) at R=11 kpc from the
Galactic center at t=5.58 Gyr).

Figure 16: The mean vertical displacement of the Galactic disc (first column)
compared to the corresponding Fourier decomposition, using the sum of the
first three terms m=0, m=1 and m=2 only (second column) and the first sixteen
terms m=0, …, m=15 (third column). Two snapshots were considered to test the
validity of the approximation both in the early stages of the simulation
(first column) and the more perturbed regimes (second row).  Figure 17: Same
as Figure 3, but considering a quiescent phase in Regime 1.  Figure 18: Same
as Figure 3, but considering Regime 2.

##  Appendix B Zoom last region

In this section, we reproduce Figures 4, 7 and 8, but zooming in the latest
and more perturbed regimes (Figure 19, 20, top and bottom, respectively). This
zooming-in can help us to identify features that otherwise would have been
missed, due to the short dynamical time-scales of Regimes 2 and 3, as
discussed in the main text.

Figure 19: Same as Figure 4 , but zooming in the most perturbed region.

Figure 20: Same as Figure 7 (top) and Figure 8 (bottom), but zooming in the
most perturbed region.

##  Appendix C Vertical spectrograms and resonances

The vertical spectrograms in Figure 10 and 11 are constructed over a temporal
baseline of 1 Gyr. Such temporal baseline has been chosen to identify
different phases, i.e. quiescent or with one (or more) interaction(s) with
Sgr, to highlight the differences in the behaviour of the Galactic disc.
However, a temporal baseline of 1 Gyr implies a resolution of 6.28 km/s/kpc on
the angular frequency (y-axis of the vertical spectrogram) (see Equation 13 in
Chequers et al. (2018)), which is quite large. Figure 21 shows, as a
complementary test, the vertical spectrogram over a baseline of 5 Gyrs, which
implies a resolution of 1.26 km/s/kpc.

Figure 22 shows vertical resonances calculated in two different ways. Red
curves are based on vertical frequencies obtained from the AGAMA package
(Vasiliev, 2018, 2019). The potentials for particles in each component in each
snapshot are approximated using a spline interpolation of radial bins combined
with fourier series (for disk components) and spherical harmonics (for
spherical components). The actions and frequencies are then calculated using
the Stackel fudge (Binney, 2012; Sanders & Binney, 2015, 2016).

Black curves show the outcome of an analytical approximation, where vertical
frequencies are written as the sum of the contribution from the halo, the
bulge and the disc
$\nu_{z}=\sqrt{\nu^{2}_{z,\rm{halo}}+\nu^{2}_{z,\rm{bulge}}+\nu^{2}_{z,\rm{disc}}}$.
Vertical frequencies from the halo and the bulge are calculated approximating
them as Hernquist profiles, evaluating Eq. 5 of D’Onghia et al. (2016) at z=0.
Following D’Onghia et al. (2013), we describe the disc density distribution as

$\rho(R,z)=\frac{M_{*}}{4\pi
z_{0}R^{2}_{s}}\rm{sech}^{2}\frac{z}{z_{0}}exp(-R/R_{s}))\quad,$ (3)

where $R_{s}$ is the scale length of the disc, $z_{0}$ is the scale height of
the disc, and $M_{*}$ is the mass of the disc. Using the above density
distribution, the vertical frequency $\nu^{2}_{z,\rm{disc}}$ can be
approximated as

$\nu^{2}_{z,\rm{disc}}=\frac{GM_{*}}{z_{0}R^{2}_{s}}\exp(-R/R_{s})\quad.$ (4)

Since the analytical approximation reflects the initial conditions, the
agreement between the black and red curves is better at the beginning of the
simulation (upper left panel) and becomes worse at the end of the simulation
(lower right panel).

Figure 21: Same as Figure 10, but showing the vertical spectrogram for the m=1
term over a time baseline of 5 Gyr (from 2 and 7 Gyr in the simulation).
Solid, dashed and dotted curves are as in Figure 10.

Figure 22: Comparison between vertical resonances from the analytical
approximation (black curves) and those obtained using vertical frequencies
from AGAMA at different times (red curves). (Top panels, from left to right:
t=1,2,3 Gyr. Bottom panels, from left to right: 4,5,6 Gyr.) Solid lines show
corotation, while the dashed and dashed-dotted curves show, respectively,
$\Omega\pm\nu_{z}$ and $\Omega\pm 2\nu_{z}$ vertical resonances.