# A Comprehensive Perturbative Formalism for Phase Mixing in Perturbed Disks.
II. Phase Spirals in an Inhomogeneous Disk Galaxy with a Non-responsive Dark
Matter Halo

Uddipan Banik Department of Astronomy, Yale University, PO. Box 208101, New
Haven, CT 06520, USA  Frank C. van den Bosch Department of Astronomy, Yale
University, PO. Box 208101, New Haven, CT 06520, USA  Martin D. Weinberg
Department of Astronomy, University of Massachusetts at Amherst, 710 N.
Pleasant St., Amherst, MA 01003

(Received XXX; Revised YYY; Accepted ZZZ)

###### Abstract

We develop a linear perturbative formalism to compute the response of an
inhomogeneous stellar disk embedded in a non-responsive dark matter (DM) halo
to various perturbations like bars, spiral arms and encounters with satellite
galaxies. Without self-gravity to reinforce it, the response of a Fourier mode
phase mixes away due to an intrinsic spread in the vertical ($\Omega_{z}$),
radial ($\Omega_{r}$) and azimuthal ($\Omega_{\phi}$) frequencies, giving rise
to local phase-space spirals. Collisional diffusion due to scattering of stars
by structures like giant molecular clouds causes super-exponential damping of
the phase-spiral amplitude. The $z-v_{z}$ phase-spiral is one-armed (two-
armed) for vertically anti-symmetric (symmetric) bending (breathing) modes.
Only transient perturbations with timescales ($\tau_{\rm P}$) comparable to
the vertical oscillation period ($\tau_{z}\sim 1/\Omega_{z}$) can trigger
vertical phase-spirals. Each $(n,l,m)$ mode of the response to impulsive
($\tau_{\rm P}<\tau=1/(n\Omega_{z}+l\Omega_{r}+m\Omega_{\phi})$) perturbations
is power law ($\sim\tau_{\rm P}/\tau$) suppressed, but that to adiabatic
($\tau_{\rm P}>\tau$) perturbations, is exponentially weak
($\sim\exp{\left[-\left(\tau_{\rm P}/\tau\right)^{\alpha}\right]}$) except
resonant ($\tau\to\infty$) modes. Slower ($\tau_{\rm P}>\tau_{z}$)
perturbations, e.g., distant encounters with satellite galaxies, induce
stronger bending modes. Sagittarius dominates the Solar neighborhood response
of the Milky Way (MW) disk to satellite encounters. Thus, if the Gaia phase-
spiral was triggered by a MW satellite, Sagittarius is the leading contender.
However, the survival of the phase-spiral against collisional damping
necessitates an impact $\sim 0.6-0.7\>{\rm Gyr}$ ago. We discuss the impact of
the detailed galactic potential on the shape of phase-spirals: phase mixing
occurs slower and thus phase-spirals are more loosely wound in the outer disk
and in presence of an ambient DM halo.

methods: analytical — Perturbation methods — Gravitational interaction —
Galaxy: disk — Galaxy: kinematics and dynamics — Galaxy stellar disks —
galaxies: interactions — Milky Way dynamics — Milky Way disk

##  1 Introduction

Disk galaxies are characterized by large-scale ordered motion and are
therefore highly responsive to perturbations. Following a time-dependent
gravitational perturbation, the actions of the disk stars are modified. This
in turn causes a perturbation in the distribution function (DF) of the disk
known as the response. Over time the response decays away as the system
‘relaxes’ towards a new quasi-equilibrium via collisionless processes that
include kinematic processes like phase mixing (loss of coherence in the
response due to different oscillation frequencies of stars) and secular/self-
gravitating/collective processes like Landau damping (loss of coherence due to
wave-particle interactions, Lynden-Bell, 1962). As pointed out by Sridhar
(1989) and Maoz (1991), phase mixing is the key ingredient of all
collisionless relaxation and re-equilibration.

The timescale of collisionless equilibration is typically longer than the
orbital periods of stars. Therefore disk galaxies usually harbour prolonged
features of incomplete equilibration following a perturbation, e.g., bars,
spiral arms, warps and other asymmetries. An intriguing example is the one-
armed phase-space spiral, or phase-spiral for short, discovered in the Gaia
DR2 data (Gaia Collaboration et al., 2018) by Antoja et al. (2018) and
discussed in more detail in subsequent studies (e.g., Bland-Hawthorn et al.,
2019; Laporte et al., 2019; Li & Widrow, 2021; Li, 2021; Gandhi et al., 2022).
Antoja et al. (2018) plotted the density of stars in the Solar neighborhood in
the $(z,v_{z})$-plane of vertical position, $z$, and vertical velocity,
$v_{z}$, and noticed a faint spiral pattern which became more pronounced when
colour-coding the $(z,v_{z})$-‘pixels’ by the median radial or azimuthal
velocities. The one-armed spiral shows 2-3 complete wraps like a snail shell,
and is interpreted as an indication of vertical phase mixing following a
perturbation that is anti-symmetric about the midplane (bending mode) and
occurred $\sim 500\>{\rm Myr}$ ago. More recently, Hunt et al. (2022) used the
more extensive Gaia DR3 data to study the distributions of stars in $z-v_{z}$
space at different locations in the MW disk. They found that unlike the one-
armed phase-spiral or bending mode at the Solar radius, the inner disk shows a
two-armed phase-spiral that corresponds to a breathing mode or symmetric
perturbation about the midplane. They inferred that while the one-armed spiral
in the Solar neighborhood might have been caused by the impact of a satellite
galaxy such as the Sagittarius dwarf, the two-armed spiral in the inner disk
could not have been induced by the same since almost all satellite impacts are
far too slow/adiabatic from the perspective of the inner disk. Rather, they
suggested that the two-armed phase-spiral might haven been triggered by a
transient spiral arm or bar.

The phase-spiral holds information about the perturbative history and
gravitational potential of the disk and can therefore serve as an essential
tool for galacto-seismology (Widrow et al., 2014; Johnston et al., 2017). For
a given potential, the winding of the spiral is an indication of the time
elapsed since the perturbation occurred with older spirals revealing more
wraps. A one-armed (two-armed) phase-spiral corresponds to a bending
(breathing) mode. Which mode dominates, in turn, depends on the time-scale of
the perturbation, with temporally shorter (longer) perturbations (e.g., a fast
or slow encounter with a satellite) predominantly triggering breathing
(bending) modes (Widrow et al., 2014; Banik et al., 2022).

In addition to depending on the nature of the perturbation, the phase-spiral
also encodes information about the oscillation frequencies of stars and thus
the detailed potential. In particular, the shape of the spiral depends on how
the vertical frequencies, $\Omega_{z}$, vary as a function of the vertical
action, $I_{z}$, which in turn depends on the underlying potential. In (Banik
et al., 2022, hereafter Paper I) we showed that the amplitude of the phase-
spiral can damp away due to lateral mixing, with a damping rate that depends
on both the spatio-temporal nature of the perturbation and the frequency
structure of the galaxy. This damping, though, only affects the response in
the coarse-grained sense, i.e., upon marginalization of the response over the
lateral degrees of freedom (the action-angle variables). Damping at the fine-
grained level requires collisional diffusion, such as that arising from the
gravitational scattering of stars against giant molecular clouds (GMCs), or
dark matter (DM) substructure (Tremaine et al., 2022).

Paper I addresses the problem of inferring the nature of the perturbation from
the amplitude and structure of the phase-spiral using a model of an infinite,
isothermal slab for the unperturbed disk. This simple, yet insightful, model
provides us with essential physical understanding of the perturbative response
of disks without the complexity of modelling a realistic, inhomogeneous disk.
However it suffers from certain glaring caveats: (i) lateral uniformity
leading to an incorrect global structure of the response in the lateral
direction, (ii) Maxwellian distribution of velocities in the lateral direction
that overpredicts lateral mixing and thereby the rate at which the amplitude
of the phase-spiral damps out, (iii) absence of a DM halo and (iv) absence of
self-gravity of the response. In this paper we relax the first three
assumptions. We consider an inhomogeneous disk characterized by a realistic DF
similar to the pseudo-isothermal DF (Binney, 2010), that properly captures the
orbital dynamics of the disk stars in 3D. In addition, we consider the effect
of an underlying DM halo which for the sake of simplicity we consider to be
non-responsive. This ambient DM halo alters the potential and thus the
frequencies of stars, which can in turn affect the shape of the phase-spiral
and its coarse-grained survival. We also consider the impact of small-scale
collisionality on the fine-grained survival of the phase-spiral. Since in this
paper we are primarily interested in the phase mixing of the disk response
that gives rise to phase-spirals, we ignore the self-gravity of the response,
which to linear order spawns coherent point mode oscillations of the disk as a
whole (for treatments of the self-gravitating response of isothermal slabs,
see Mathur, 1990; Weinberg, 1991) and somewhat enhances the amplitude of
phase-spirals.

This paper is organized as follows. Section 2 describes the standard linear
perturbation theory for collisionless systems and its application to a
realistic disk galaxy embedded in a DM halo that is exposed to a general
perturbation. Sections 3 and 4 are concerned with computing the disk response
for different perturber models. In Section 3 we compute the disk response and
phase-spirals induced by bars and spiral arms. We also discuss the impact of
collisional diffusion on the fine-grained survivability of the phase-spiral.
In Section 4 we compute the response to encounters with satellite galaxies.
Section 5 describes how phase-spirals can be used to constrain the galactic
potential. We summarize our findings in Section 6.

##  2 Linear perturbation theory for galaxies

###  2.1 Linear perturbative formalism

A galaxy, to very good approximation, is devoid of star-star collisions.
However, there are other potential sources of collisions such as scatterings
due to gravitational interactions of stars with giant molecular clouds (GMCs)
or DM substructure. The dynamics of stars in such a system is governed by the
Boltzmann equation:

$\displaystyle\frac{\partial f}{\partial t}+[f,H]=C[f],$ (1)

where $f$ denotes the DF, $H$ denotes the Hamiltonian, square brackets denote
the Poisson bracket, and $C[f]$ denotes the collision operator due to small-
scale fluctuations, which can be approximated by a Fokker-Planck operator (see
Appendix A of Tremaine et al., 2022):

$\displaystyle
C[f]=\frac{1}{2}\frac{\partial}{\partial\xi_{i}}\left(D_{ij}\frac{\partial
f}{\partial\xi_{j}}\right),$ (2)

where $\bm{\xi}=\left({\bf q},{\bf p}\right)$ with ${\bf q}$ and ${\bf p}$
denoting the canonically conjugate position and momentum variables, and
$D_{ij}$ denotes the diffusion coefficient tensor.

Let the unperturbed steady state Hamiltonian of the galaxy be $H_{0}$ and the
corresponding DF be given by $f_{0}$, which satisfies the unperturbed Fokker-
Planck equation (FPE),

$\displaystyle[f_{0},H_{0}]=C[f_{0}].$ (3)

In presence of a small time-dependent perturbation in the potential,
$\Phi_{\rm P}(t)$, the perturbed Hamiltonian can be written as

$\displaystyle H=H_{0}+\Phi_{\rm P}(t)+\Phi_{1}(t),$ (4)

where $\Phi_{1}$ is the gravitational potential related to the response
density, $\rho_{1}=\int f_{1}{\rm d}^{3}{\bf v}$, via the Poisson equation,

$\displaystyle\nabla^{2}\Phi_{1}=4\pi G\rho_{1}.$ (5)

The perturbed DF can be written as

$\displaystyle f=f_{0}+f_{1},$ (6)

where $f_{1}$ is the linear order perturbation in the DF. In the weak
perturbation limit where linear perturbation theory holds, the time-evolution
of $f_{1}$ is dictated by the following linearized form of the FPE:

$\displaystyle\frac{\partial f_{1}}{\partial t}+[f_{1},H_{0}]+[f_{0},\Phi_{\rm
P}]+[f_{0},\Phi_{1}]=C[f_{1}].$ (7)

Throughout this paper we neglect the self-gravity of the disk response, which
implies that we set the polarization term, $[f_{0},\Phi_{1}]=0$. The
implications of including self-gravity are discussed in Paper I.

###  2.2 Response of a Galactic Disk to a realistic perturbation

The dynamics of a realistic disk galaxy like the Milky Way (MW) is quasi-
periodic, i.e., can be characterized by oscillations in the azimuthal, radial
and vertical directions. In close proximity to the mid-plane and under radial
epicyclic approximation, the Hamilton-Jacobi equation becomes separable,
implying that all stars confined within a few vertical scale heights from the
mid-plane of the disk are on regular, quasi-periodic orbits that are
characterized by a radial action, $I_{R}$, an azimuthal action $I_{\phi}$, and
a vertical action $I_{z}$. Hence, the motion of each star is characterized by
three frequencies:

$\displaystyle\Omega_{R}=\frac{\partial H_{0}}{\partial
I_{R}}\,,\;\;\;\Omega_{\phi}=\frac{\partial H_{0}}{\partial
I_{\phi}}\,,\;\;\;\Omega_{z}=\frac{\partial H_{0}}{\partial I_{z}}\,.$ (8)

This quasi-periodic nature of the orbits near the mid-plane is approximately
preserved even in the presence of a (non-triaxial) DM halo since this
preserves the axi-symmetry of the potential. Typically, as discussed in
section 5, the presence of a halo increases the oscillation frequencies of the
disk stars.

In terms of these canonical conjugate action-angle variables, using equation
(8), the linearized form of the FPE given in Equation (7) becomes

$\displaystyle\frac{\partial f_{1}}{\partial t}+\Omega_{z}\frac{\partial
f_{1}}{\partial w_{z}}+\Omega_{R}\frac{\partial f_{1}}{\partial
w_{R}}+\Omega_{\phi}\frac{\partial f_{1}}{\partial
w_{\phi}}-\frac{\partial\Phi_{\rm P}}{\partial w_{z}}\frac{\partial
f_{0}}{\partial I_{z}}-\frac{\partial\Phi_{\rm P}}{\partial
w_{R}}\frac{\partial f_{0}}{\partial I_{R}}-\frac{\partial\Phi_{\rm
P}}{\partial w_{\phi}}\frac{\partial f_{0}}{\partial I_{\phi}}$
$\displaystyle=D_{z}\frac{\partial}{\partial I_{z}}\left(I_{z}\frac{\partial
f_{1}}{\partial I_{z}}\right)$
$\displaystyle+\frac{D_{z}}{4I_{z}}\frac{\partial^{2}f_{1}}{\partial
w^{2}_{z}}+\frac{D_{R}}{4I_{R}}\frac{\partial^{2}f_{1}}{\partial w^{2}_{R}}$
(9)

Here we have performed certain simplifications of the Fokker-Planck operator.
Firstly, following Binney & Lacey (1988), we have assumed that
$D_{zz}=D_{z}I_{z}$ and $D_{RR}=D_{R}I_{R}$, since this preserves the pseudo-
isothermal form of the unperturbed DF (equation [15]) of the disk. Secondly,
the $I_{R}$ diffusion of the response $f_{1}$ is negligible. This is because
the frequencies do not depend on $I_{R}$ under the radial epicyclic
approximation (and only mildly depend on $I_{R}$ without it) and therefore the
response does not develop $I_{R}$ gradients. Thirdly, following Binney & Lacey
(1988), we have neglected diffusion in $I_{\phi}$ and $w_{\phi}$ since the
terms involving $D_{\phi\phi}$, $D_{r\phi}$ and $D_{\phi z}$ are smaller than
the $I_{z}$ and $I_{R}$ diffusion terms by factors of at least
$\sigma_{R}/v_{c}$ or $\sigma_{z}/v_{c}$, which are typically much smaller
than unity ($\sigma_{R}$ and $\sigma_{z}$ are radial and vertical velocity
dispersions respectively, and $v_{c}$ is the circular velocity along $\phi$).
We have retained the $w_{z}$ and $w_{R}$ diffusion terms for the sake of
completeness, but as we point out later, the diffusion in angles typically
occurs over much longer timescales than that in actions and hence is
comparatively less important.

Since the stars move along quasi-periodic orbits characterized by actions and
angles, we can expand the perturbations, $\Phi_{\rm P}$ and $f_{1}$, as
discrete Fourier series in the angles as follows

$\displaystyle\Phi_{\rm P}\left({\bf w},{\bf I},t\right)$
$\displaystyle=\sum_{n=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}\exp{\left[i(nw_{z}+lw_{R}+mw_{\phi})\right]}\,\Phi_{nlm}\left({\bf
I},t\right),$ $\displaystyle f_{1}\left({\bf w},{\bf I},t\right)$
$\displaystyle=\sum_{n=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}\exp{\left[i(nw_{z}+lw_{R}+mw_{\phi})\right]}\,f_{1nlm}({\bf
I},t),$ (10)

where ${\bf w}=(w_{z},w_{R},w_{\phi})$ and ${\bf I}=(I_{z},I_{R},I_{\phi})$.
Substituting these Fourier expansions in equation (9) yields the following
differential equation for the evolution of $f_{1nlm}$:

$\displaystyle\frac{\partial f_{1nlm}}{\partial
t}+i(n\Omega_{z}+l\Omega_{R}+m\Omega_{\phi})f_{1nlm}$
$\displaystyle=i\left(n\frac{\partial f_{0}}{\partial I_{z}}+l\frac{\partial
f_{0}}{\partial I_{R}}+m\frac{\partial f_{0}}{\partial
I_{\phi}}\right)\Phi_{nlm}$ $\displaystyle+D_{z}\frac{\partial}{\partial
I_{z}}\left(I_{z}\frac{\partial f_{1nlm}}{\partial
I_{z}}\right)-\left[\frac{n^{2}D_{z}}{4I_{z}}+\frac{l^{2}D_{R}}{4I_{R}}\right]f_{1nlm}.$
(11)

This can be solved using the Green’s function technique, with the initial
condition, $f_{1nlm}(t_{\rm i})=0$, to yield the following closed integral
form for $f_{1nlm}$:

$\displaystyle f_{1nlm}({\bf I},t)$ $\displaystyle=i\left(n\frac{\partial
f_{0}}{\partial I_{z}}+l\frac{\partial f_{0}}{\partial I_{R}}+m\frac{\partial
f_{0}}{\partial I_{\phi}}\right){\cal I}_{nlm}({\bf I},t).$ (12)

Here, for $D_{z}\ll\sigma^{2}_{z}$ ($\sigma_{z}$ is the vertical velocity
dispersion), which is typically the case, ${\cal I}_{nlm}({\bf I},t)$ can be
approximately expressed as

$\displaystyle{\cal I}_{nlm}({\bf I},t)$ $\displaystyle\approx\int_{t_{\rm
i}}^{t}{\rm d}\tau\,{\cal G}_{nlm}({\bf I},t-\tau)\,\Phi_{nlm}({\bf I},\tau).$
(13)

Here ${\cal G}_{nlm}(t-\tau)$ is the Green’s function (see Appendix A of
Tremaine et al. (2022) for derivation), given by

$\displaystyle\mathcal{G}_{nlm}({\bf I},t-\tau)$
$\displaystyle\approx\exp{\left[-i(n\Omega_{z}+l\Omega_{R}+m\Omega_{\phi})(t-\tau)\right]}$
$\displaystyle\times\exp{\left[-\left(\frac{n^{2}D_{z}}{4I_{z}}+\frac{l^{2}D_{R}}{4I_{R}}\right)\left(t-\tau\right)\right]}\,\exp{\left[-\frac{{\left(n\Omega_{z1}\right)}^{2}D_{z}I_{z}}{3}{\left(t-\tau\right)}^{3}\right]},$
(14)

where $\Omega_{z1}=\partial\Omega_{z}/\partial I_{z}$. The sinusoidal factor
represents the oscillations of stars at their natural frequencies which vary
with actions, leading to the formation of phase-spirals (see section 3.1.1 for
details). The first exponential damping factor indicates the damping of the
response due to diffusion in angles while the second damping factor manifests
the damping of the $I_{z}$ gradients of the response by diffusion in $I_{z}$.
As discussed in section 3.1.2, the diffusion in actions is much more efficient
than that in angles.

Each $(n,l,m)$ Fourier coefficient of the response acts as a forced damped
oscillator with three different natural frequencies, $n\Omega_{z}$,
$l\Omega_{R}$ and $m\Omega_{\phi}$, which is being driven by an external time-
dependent perturber potential, $\Phi_{nlm}$, and damped due to collisional
diffusion. A similar expression, albeit without allowing for collisionality,
for the DF perturbation has been derived by Carlberg & Sellwood (1985) in the
context of spiral arm induced perturbations and radial migrations in the
galactic disk, and by other previous studies (e.g., Lynden-Bell & Kalnajs,
1972; Tremaine & Weinberg, 1984; Carlberg & Sellwood, 1985; Weinberg, 1989,
1991, 2004; Kaur & Sridhar, 2018; Banik & van den Bosch, 2021a; Kaur & Stone,
2022) in the context of dynamical friction in spherical systems. To obtain the
final expression for $f_{1nlm}$, we need to specify the DF $f_{0}$ of the
unperturbed galaxy, as well as the spatio-temporal behavior of the perturber
potential, $\Phi_{\rm P}$, which is addressed below.

###  2.3 The unperturbed galaxy

Under the radial epicyclic approximation (small $I_{R}$), the unperturbed DF,
$f_{0}$, for a rotating MW-like disk galaxy can be well approximated as a
pseudo-isothermal DF, i.e., written as a nearly isothermal separable function
of the azimuthal, radial and vertical actions. Following Binney (2010), we
write

$\displaystyle
f_{0}\approx\frac{\sqrt{2}}{\pi^{3/2}\,\sigma_{z}h_{z}}{\left(\frac{\Omega_{\phi}\Sigma}{\kappa\,\sigma^{2}_{R}}\right)}_{R_{\rm
c}}\,\exp{\left[-\frac{\kappa
I_{R}}{\sigma^{2}_{R}}\right]}\,\exp{\left[-\frac{E_{z}(I_{z})}{\sigma^{2}_{z}}\right]}\,\Theta(L_{z})\,,$
(15)

The vertical structure of this disk is isothermal, while the radial profile is
pseudo-isothermal. Here
$\Sigma=\Sigma(R)=\int_{-\infty}^{\infty}dz\,\rho(R,z)$ is the surface density
of the disk, $L_{z}$ is the $z$-component of the angular momentum, which is
equal to $I_{\phi}$, $R_{\rm c}=R_{\rm c}(L_{z})$ is the guiding radius,
$\Omega_{\phi}$ is the circular frequency, and $\kappa=\kappa(R_{\rm
c})=\lim_{I_{R}\to 0}{\Omega_{R}}$ is the radial epicyclic frequency (Binney &
Tremaine, 1987). $\Theta(x)$ is the Heaviside step function. Thus we assume
that the entire galaxy is composed of prograde stars with $L_{z}>0$.

The density profile, $\rho(R,z)$, of the disk corresponding to the above DF is
the product of a radially exponential profile with scale radius $h_{r}$ and a
vertically isothermal ($\operatorname{sech}^{2}$) profile with scale height
$h_{z}$ (equation [A3]). As shown by Smith et al. (2015), this density profile
is accurately approximated by a sum of three Miyamoto & Nagai (1975)
disks111the 3MN profile as implemented in the Gala Python package (Price-
Whelan, 2017; Price-Whelan et al., 2020)., which has a simple, analytical form
for the associated potential. Throughout, we therefore use this 3MN
approximation for our disk since this drastically simplifies the computation
of orbital frequencies. The disk is assumed to be embedded in an extended DM
halo characterized by a spherical NFW (Navarro et al., 1997) density profile,
with virial mass $M_{\rm vir}$, concentration $c$, scale radius $r_{s}$ and
the corresponding potential $\Phi_{\rm h}$ given by equation (A6). Throughout,
for the purpose of computing the disk response, we assume typical MW like
parameters for the various quantities, i.e., $\>{R_{\odot}}=8$ kpc, disk mass
$M_{\rm d}=5\times 10^{10}\>{\rm M_{\odot}}$, $h_{R}=2.2\>{\rm kpc}$,
$\sigma_{R}(\>{R_{\odot}})=\sigma_{R,\odot}=35\,{\rm km}/{\rm s}$,
$h_{z}=0.4\>{\rm kpc}$ and
$\sigma_{z}(\>{R_{\odot}})=\sigma_{z,\odot}=\sqrt{2\pi
Gh_{z}\Sigma(\>{R_{\odot}})}=23\,{\rm km}/{\rm s}$ (McMillan, 2011; Bovy &
Rix, 2013). For the NFW DM halo, we adopt $M_{\rm vir}=9.78\times
10^{11}\>{\rm M_{\odot}}$, $r_{s}=16$ kpc, and $c=15.3$ (Bovy, 2015).

The combined potential experienced by the disk stars is simply the sum of disk
and halo potentials, i.e.,

$\displaystyle\Phi_{0}(R,z)=\Phi_{\rm d}(R,z)+\Phi_{\rm h}(R,z).$ (16)

The total energy of a disk star under the radial epicyclic approximation is
$E=L^{2}_{z}/2R^{2}_{c}+\Phi_{0}(R_{\rm c},0)+\kappa I_{R}+E_{z}$, where the
vertical part of the energy is given by $E_{z}=v^{2}_{z}/2+\Phi_{z}(R_{\rm
c},z)$, with $R_{\rm c}(L_{z})$ the guiding radius given by $L^{2}_{z}/R_{\rm
c}^{3}=\partial\Phi_{0}/\partial R|_{R=R_{\rm c}}$. The vertical potential,
$\Phi_{z}(R_{\rm c},z)$, is given by

$\displaystyle\Phi_{z}(R_{\rm c},z)=\Phi_{0}(R_{\rm c},z)-\Phi_{0}(R_{\rm
c},0).$ (17)

The vertical action, $I_{z}$, can be obtained from $E_{z}$ as follows

$\displaystyle I_{z}=\frac{1}{2\pi}\oint
v_{z}\,dz=\frac{2}{\pi}\int_{0}^{z_{\rm max}}\sqrt{2[E_{z}-\Phi_{z}(R_{\rm
c},z)]}\,dz,$ (18)

where $\Phi_{z}(R_{\rm c},z_{\rm max})=E_{z}$. This implicit equation can be
inverted to obtain $E_{z}(R_{\rm c},I_{z})$. The time period of vertical
oscillation can then be obtained using

$\displaystyle T_{z}(R_{\rm c},I_{z})=\oint\frac{dz}{v_{z}}=4\int_{0}^{z_{\rm
max}}\frac{dz}{\sqrt{2\left[E_{z}(R_{\rm c},I_{z})-\Phi_{z}(R_{\rm
c},z)\right]}},$ (19)

which yields the vertical frequency, $\Omega_{z}(R_{\rm
c},I_{z})=2\pi/T_{z}(R_{\rm c},I_{z})$.

Substituting the expression for $f_{0}$ given by Equation (15) in Equation
(12), we obtain the following integral form for $f_{1nlm}$,

$\displaystyle f_{1nlm}({\bf I},t)$
$\displaystyle\approx-\frac{2i}{\pi\sigma^{2}_{R}}\,\frac{1}{\sqrt{2\pi}h_{z}\sigma_{z}}\exp{\left[-\frac{\kappa
I_{R}}{\sigma^{2}_{R}}\right]}\exp{\left[-\frac{E_{z}(I_{z})}{\sigma^{2}_{z}}\right]}$
$\displaystyle\times\left[\left\\{\left(\frac{n\Omega_{z}}{\sigma^{2}_{z}}+\frac{l\kappa}{\sigma^{2}_{R}}\right){\left(\frac{\Omega_{\phi}\Sigma}{\kappa}\right)}-m\frac{{\rm
d}}{{\rm
d}L_{z}}\left(\frac{\Omega_{\phi}\Sigma}{\kappa}\right)\right\\}\Theta(L_{z})-m\frac{\Omega_{\phi}\Sigma}{\kappa}\delta(L_{z})\right]{\cal
I}_{nlm}({\bf I},t).$ (20)

As we shall see, the first order disk response expressed above phase mixes
away and gives rise to phase-spirals due to oscillations of stars with
different frequencies except when they are resonant with the frequency of the
perturber. However this ‘direct’ response of the disk does not include certain
effects. First of all, we ignore the self-gravity of the response. As
discussed in Paper I, to linear order self-gravity gives rise to point mode
oscillations of the disk that are decoupled from the phase mixing component of
the response which is what we are interested in. Secondly, for the sake of
simplicity, we consider the ambient DM halo to be non-responsive and therefore
ignore the indirect effect of the halo response on disk oscillations. We leave
the inclusion of these two effects in the computation of the disk response for
future work.

The spatio-temporal nature of the perturbing potential dictates the disk
response. In this paper we explore two different types of perturbation to
which realistic disc galaxies can be exposed, and which are thus of general
astrophysical interest. The first is an in-plane spiral/bar perturbation with
a vertical structure, either formed as a consequence of secular evolution, or
triggered by an external perturbation. We consider both short-lived
(transient) and persistent spirals. The second type of perturbation that we
consider is that due to an encounter with a massive object, e.g., a satellite
galaxy or DM subhalo.

##  3 Disk response to spiral arms and bars

We model the potential of a spiral arm perturbation as one with a vertical
profile and a sinusoidal variation along radial and azimuthal directions,

$\displaystyle\Phi_{\rm P}(R,\phi,z)$ $\displaystyle=-\frac{2\pi G\Sigma_{\rm
P}}{k_{R}}\,\left[\alpha\,{\cal M}_{\rm o}(t)\,{\cal F}_{\rm o}(z)+{\cal
M}_{\rm e}(t)\,{\cal F}_{\rm
e}(z)\right]\sum_{m_{\phi}=0,2}\sin{\left[k_{R}R+m_{\phi}\left(\phi-\Omega_{\rm
P}t\right)\right]}\,.$ (21)

Here $\Omega_{\rm P}$ is the pattern speed and $k_{R}$ is the horizontal wave
number of the spiral perturbation. The long wavelength limit, $k_{R}\to 0$,
corresponds to a bar. We consider the in-plane part of $\Phi_{\rm P}$ to be a
combination of an axisymmetric ($m_{\phi}=0$) and a 2-armed spiral mode
($m_{\phi}=2$), and the vertical part to be a combination of anti-
symmetric/odd and symmetric/even perturbations respectively denoted by ${\cal
F}_{\rm o}$ and ${\cal F}_{\rm e}$, that are modulated by growth functions,
${\cal M}_{\rm o}(t)$ and ${\cal M}_{\rm e}(t)$, capturing the growth and/or
decay of the spiral strength over time. The ratio of the maximum strengths of
the anti-symmetric and symmetric parts of the perturbation is $\alpha$. We
consider the following two functional forms for ${\cal M}_{j}(t)$ (where the
subscript $j={\rm o}$ or ${\rm e}$):

$\displaystyle{\cal M}_{j}(t)$
$\displaystyle=\begin{cases}\frac{1}{\sqrt{\pi}}\exp{\left[-\omega^{2}_{j}t^{2}\right]},&\text{Transient
spiral/bar}\\\
\exp{\left[\gamma_{j}t\right]}+\left(1-\exp{\left[\gamma_{j}t\right]}\right)\,\Theta(t),&\text{Persistent
spiral/bar.}\end{cases}$ (22)

The first option describes a transient spiral/bar that grows and decays like a
Gaussian pulse with a characteristic life-time $\tau_{{\rm P}j}\sim
1/\omega_{j}$ (Banik et al., 2022). The second form describes a persistent
spiral perturbation that grows exponentially on a timescale $\tau_{{\rm
G}j}\sim 1/\gamma_{j}$ and then saturates to a constant amplitude. We shall
see shortly that these two kinds of spiral perturbations perturb the disk in
very different ways.

The vertical part of the perturbation consists of an anti-symmetric function,
${\cal F}_{\rm o}(z)$, and a symmetric function, ${\cal F}_{\rm e}(z)$, which,
for the sake of simplicity, we take to be the following trigonometric
functions:

$\displaystyle{\cal F}_{\rm o}(z)$ $\displaystyle=\sin{\left(k_{z}^{({\rm
o})}z\right)},$ $\displaystyle{\cal F}_{\rm e}(z)$
$\displaystyle=\cos{\left(k_{z}^{({\rm e})}z\right)}.$ (23)

Here $k^{({\rm o})}_{z}$ and $k^{({\rm e})}_{z}$ denote the vertical wave-
numbers of the anti-symmetric and symmetric perturbations, respectively. Since
the above functions form a complete Fourier basis in $z$, any (vertical)
perturber profile can be expressed as a linear superposition of ${\cal F}_{\rm
o}$ and ${\cal F}_{\rm e}$. The disk response involves the Fourier
coefficients of the perturbing potential, $\Phi_{nlm}$, which can be obtained
by taking the Fourier transform of $\Phi_{\rm P}$ given in Equation (21) with
respect to the angles, $w_{R}$, $w_{\phi}$ and $w_{z}$, as detailed in
Appendix B.

Figure 1: MW disk response to transient bars/2-armed spirals with Gaussian
temporal modulation in absence of collisional diffusion: Left panel shows the
steady state ($t\to\infty$) amplitude of the disk response, $f_{1,nlm}/f_{0}$,
in the Solar neighborhood, computed using equations (24) and (27) in presence
of an ambient DM halo, as a function of the pulse frequency, $\omega_{j}$,
where the subscript $j={\rm o}$ and ${\rm e}$ for vertically anti-symmetric
(odd $n$) and symmetric (even $n$) perturbations. Solid (dashed) lines
indicate the $n=1$ bending ($n=2$ breathing) modes and different colors denote
$(l,m)=(0,-2)$, $(0,0)$ and $(0,2)$ respectively. We consider
$I_{z}=I_{z,\odot}\equiv h_{z}\sigma_{z,\odot}$ and marginalize the response
over $I_{R}$. Note that the response peaks at intermediate values of
$\omega_{j}$, which is different for different modes, and is suppressed like a
power law in the impulsive (large $\omega_{j}$) limit and super-exponentially
in the adiabatic (small $\omega_{j}$) limit. Right panel shows the breathing-
to-bending ratio, $f_{1,200}/f_{1,100}$, as a function of $\omega_{\rm e}$ and
$\omega_{\rm o}$, the pulse frequencies of the bending and breathing mode
perturbations respectively. The dashed, solid, dot-dashed and dotted contours
correspond to breathing-to-bending ratios of $0.1,1,5$ and $10$ respectively.
The breathing-to-bending ratio rises and falls with increasing $\omega_{\rm
e}$ at fixed $\omega_{\rm o}$, while the reverse occurs with increasing
$\omega_{\rm o}$ at fixed $\omega_{\rm e}$, leading to a saddle point at
$(\omega_{\rm e},\omega_{\rm o})\approx(9,7)$.

###  3.1 Computing the disk response

The expression for the disk response to bars or spiral arms can be obtained by
substituting the Fourier coefficient of the perturber potential given in
Equation (B5) in Equation (20) and performing the $\tau$ integration with the
initial time, $t_{\rm i}\to-\infty$. This yields the modal response,
$f_{1nlm}$ (Equation [20]), with ${\cal I}_{nlm}({\bf I},t)$ given by

$\displaystyle{\cal I}_{nlm}({\bf I},t)$ $\displaystyle=\alpha\,\Psi^{({\rm
o})}_{nlm}({\bf I})\,{\cal P}^{({\rm o})}_{nlm}({\bf I},t)+\Psi^{({\rm
e})}_{nlm}({\bf I})\,{\cal P}^{({\rm e})}_{nlm}({\bf I},t),$ (24)

where $\Psi^{({\rm o})}_{nlm}$ and $\Psi^{({\rm e})}_{nlm}$ respectively
denote the time-independent parts of the odd and even terms in the expression
for $\Phi_{nlm}$, and

$\displaystyle{\cal P}^{(j)}_{nlm}({\bf I},t)$
$\displaystyle=\exp{\left[-i\,m\Omega_{\rm P}\,t\right]}$
$\displaystyle\times\int_{0}^{\infty}{\rm d}\tau\exp{\left[-i\,\Omega_{\rm
res}\,\tau\right]}\,\exp{\left[-\left(\frac{n^{2}D_{z}}{4I_{z}}+\frac{l^{2}D_{R}}{4I_{R}}\right)\tau\right]}\,\exp{\left[-\frac{{\left(n\Omega_{z1}\right)}^{2}D_{z}I_{z}}{3}{\tau}^{3}\right]}\,{\cal
M}_{j}(t-\tau),$ (25)

which characterizes the temporal evolution of the response. Here the subscript
$j={\rm o}$ or ${\rm e}$, and the resonance frequency, $\Omega_{\rm res}$, is
given by

$\displaystyle\Omega_{\rm res}=n\Omega_{z}+l\kappa+m(\Omega_{\phi}-\Omega_{\rm
P}).$ (26)

####  3.1.1 Collisionless limit

First we examine the response in the limit of zero diffusion, i.e., $D_{z}=0$,
where each star acts as a forced oscillator.

##### Transient spirals and bars

First we consider the case of transient spiral arm or bar perturbations that
grow and decay in strength over time, i.e., the temporal modulation ${\cal
M}_{j}(t)$ is given by the first of equations (22). In this case,

$\displaystyle{\cal P}_{nlm}^{(j)}({\bf I},t)$
$\displaystyle=\frac{1}{2\,\omega_{j}}\,\exp{\left[-\frac{\Omega^{2}_{\rm
res}}{4\omega^{2}_{j}}\right]}\left[1+\operatorname{erf}{\left(\omega_{j}t-i\frac{\Omega_{\rm
res}}{2\omega_{j}}\right)}\right]\exp{\left[-i(n\Omega_{z}+l\kappa+m\Omega_{\phi})t\right]}$
$\displaystyle\xrightarrow{t\to\infty}\frac{1}{\omega_{j}}\,\exp{\left[-\frac{\Omega^{2}_{\rm
res}}{4\omega^{2}_{j}}\right]}\,\exp{\left[-i(n\Omega_{z}+l\kappa+m\Omega_{\phi})t\right]}.$
(27)

The error function describes the growth and transient oscillations of the
response amplitude; over time the transients die away, and in the limit
$t\to\infty$ the response saturates to a constant amplitude (in the absence of
collisional diffusion).

The left-hand panel of Fig. 1 plots the amplitude of the steady state disk
response to transient spiral/bar perturbations, relative to the unperturbed
DF, as a function of the modulation/pulse frequency, $\omega_{j}$ ($j={\rm o}$
and ${\rm e}$ for bending and breathing modes respectively), for different
modes indicated in different colors. Solid and dashed lines correspond to the
$n=1$ bending modes and the $n=2$ breathing modes, respectively. We adopt
$\Sigma_{\rm P}=5.5\>{\rm M_{\odot}}\>{\rm pc}^{-2}$, $\Omega_{\rm P}=12\>{\rm
km}\,{\rm s}^{-1}\>{\rm kpc}^{-1}$, $k^{(o)}_{z}=k^{(e)}_{z}=1\>{\rm kpc}$,
$k_{R}=10\>{\rm kpc}$, and $I_{z}=I_{z,\odot}\equiv
h_{z}\sigma_{z,\odot}=9.2\>{\rm km}\,{\rm s}^{-1}$, and marginalize the
response over $I_{R}$. We set $\alpha=1$, implying equal maximum strengths for
the bending and breathing modes. As evident from this figure, and also from
equation (27), the long-term strength of the disk response (after the initial
transients have died out like $e^{-\omega^{2}_{j}t^{2}}$) scales as $\sim
1/\omega_{j}$ in the impulsive (large $\omega_{j}$) limit, but is super-
exponentially suppressed ($\sim\exp\left[-\Omega^{2}_{\rm
res}/4\omega^{2}_{j}\right]$) in the adiabatic (small $\omega_{j}$) limit away
from resonances, i.e., for $\Omega_{\rm res}\neq 0$. The adiabatic suppression
scales differently with $\omega_{j}$ for other functional forms of ${\cal
M}_{j}(t)$, e.g., for ${\cal M}_{j}(t)=1/\sqrt{1+\omega^{2}_{j}t^{2}}$ the
response strength is exponentially suppressed $(\sim\exp[-\Omega_{\rm
res}/\omega_{j}])$. The response of resonant modes ($\Omega_{\rm res}=0$)
however does not undergo adiabatic suppression and scales as $\sim
1/\omega_{j}$ throughout, becoming non-linear in the adiabatic regime. Since
there are many resonance modes, the cumulative response in the adiabatic limit
of all modes combined is only suppressed as a power-law, rather than an
exponential, in $\omega_{j}$ (Weinberg, 1994a, b).

The sinusoidal factor,
$\exp{\left[-i(n\Omega_{z}+l\kappa+m\Omega_{\phi})t\right]}$, in ${\cal
P}_{nlm}^{(j)}$ describes the oscillations of stars with three different
frequencies, $\Omega_{z}$, $\kappa$ and $\Omega_{\phi}$, along the vertical,
radial and azimuthal directions, respectively. Due to the dependence of these
frequencies on the actions, that of $\Omega_{z}$ on $I_{z}$ and of $\kappa$
and $\Omega_{\phi}$ on $I_{\phi}=L_{z}$, the response integrated over actions
eventually phase mixes away. This manifests as phase-spirals in the
$I_{z}\cos{w_{z}}-I_{z}\sin{w_{z}}$ and
$I_{\phi}\cos{\phi}-I_{\phi}\sin{\phi}$ phase-spaces, which are proxies for
the $z-v_{z}$ and $\phi-\dot{\phi}$ phase-spaces, respectively. As is evident
from equation (25), ${\cal P}_{nlm}^{(j)}\sim\exp{\left[-im\Omega_{\rm
P}t\right]}$ in the adiabatic limit ($\omega_{j}\to 0$); hence, in this limit
the sinusoidal factor,
$\exp{\left[-i(n\Omega_{z}+l\kappa+m\Omega_{\phi})t\right]}$ is absent from
the response, which implies that phase-spirals only occur for sufficiently
impulsive perturbations. As shown in Paper I, $n=1$ bending modes involve a
dipolar perturbation in the vertical phase-space
($I_{z}\cos{w_{z}}-I_{z}\sin{w_{z}}$) distribution immediately after the
perturbing pulse reaches its maximum strength. This dipolar distortion is
subsequently wound up into a one-armed phase-spiral since $\Omega_{z}$ is a
function of $I_{z}$. Breathing modes, on the other hand, involve an initial
quadrupolar perturbation in the phase-space distribution which is subsequently
wrapped up into a two-armed phase-spiral. Since $\Omega_{z}$, $\Omega_{\phi}$
and $\Omega_{R}$ all depend on $L_{z}$, the amplitude of the
$I_{z}\cos{w_{z}}-I_{z}\sin{w_{z}}$ phase-spiral damps out over time due to
mixing between stars with different $L_{z}$. The modal response, $f_{1nlm}$,
when marginalized in a narrow bin of size $\Delta L_{z}$ around $L_{z}$, damps
out as follows:

$\displaystyle\left<f_{1nlm}\right>({\bf I},t)$ $\displaystyle=\frac{1}{\Delta
L_{z}}\int_{L_{z}-\Delta L_{z}/2}^{L_{z}+\Delta L_{z}/2}{\rm
d}L_{z}\,f_{1nlm}({\bf I},t)$
$\displaystyle\approx\dfrac{\;\;\;\sin{\left[\left(\dfrac{\partial}{\partial
L_{z}}\left(n\Omega_{z}+l\kappa+m\Omega_{\phi}\right)\right)\dfrac{\Delta
L_{z}}{2}t\right]}\;\;\;}{\;\;\;\;\;\left(\dfrac{\partial}{\partial
L_{z}}\left(n\Omega_{z}+l\kappa+m\Omega_{\phi}\right)\right)\dfrac{\Delta
L_{z}}{2}t}\,f_{1nlm}({\bf I},t).$ (28)

Since the frequencies vary with $L_{z}$, marginalizing over $L_{z}$ mixes
phase-spirals that differ slightly in phases, giving way to a $\sim 1/t$
damping accompanied by a beat-like modulation with a characteristic lateral
mixing timescale,

$\displaystyle\tau_{\rm D}^{(\rm{LM})}$
$\displaystyle=\dfrac{1}{\left(\dfrac{\partial}{\partial
L_{z}}\left(n\Omega_{z}+l\kappa+m\Omega_{\phi}\right)\right)\dfrac{\Delta
L_{z}}{2}}.$ (29)

This explains why the density-contrast of the Gaia phase-spiral is enhanced
upon color-coding by $v_{\phi}$ or, equivalently, $L_{z}$ (Antoja et al.,
2018; Bland-Hawthorn et al., 2019). Radial phase mixing is also present, but
is typically much weaker because none of the frequencies depend on $I_{R}$
under the radial epicyclic approximation and only mildly depend on $I_{R}$
without it. Hence, due to ordered motion, the phase-spiral amplitude in a
realistic disk galaxy damps out at a much slower rate, as $\sim 1/t$ (in
absence of collisional diffusion), than the lateral mixing damping in the
isothermal slab case considered in Paper I, which arises from the
unconstrained lateral velocities of the stars and exhibits a Gaussian temporal
behavior.

It is worth emphasizing that not all frequencies undergo phase mixing. In fact
the resonant frequencies, for which

$\displaystyle\Omega_{\rm res}=n\Omega_{z}+l\kappa+m(\Omega_{\phi}-\Omega_{\rm
P})=0,$ (30)

do not phase mix away. Hence, parts of the phase-space closer to a resonance
take longer to phase-mix away. Moreover, as manifest from the adiabatic
suppression factor, $\exp[-\Omega^{2}_{\rm res}/4\omega^{2}_{j}]$, the near-
resonant modes with $\Omega_{\rm res}\ll 2\omega_{j}$ have much larger
amplitude than those with $\Omega_{\rm res}\gg 2\omega_{j}$ that are far from
resonance. Therefore the long-term disk response consists of stars in (near)
resonance with the perturbing bar or spiral arm. Most of the strong resonances
are confined to the disk-plane, including the co-rotation resonance
$(n,l,m)=(0,0,m)$, the Lindblad resonances $(0,\pm 1,\pm 2)$, the
ultraharmonic resonances $(0,\pm 1,\pm 4)$, and so on. For thin disks with
$h_{z}\ll h_{R}$ , the vertical degrees of freedom are generally not in
resonance with the radial or azimuthal ones since $\Omega_{z}$ is much larger
than $\Omega_{\phi}$ or $\kappa$. Hence the vertical oscillation modes ($n\neq
0$) such as the $n=1$ bending or $n=2$ breathing modes undergo phase mixing
and give rise to phase-spirals. However, if the disk has significant
thickness, then the vertical degrees of freedom can be in resonance with the
horizontal ones, e.g., banana orbits ($\Omega_{z}=2\Omega_{r}$) in barred
disks.

The excitability of the bending and breathing modes is dictated by the
perturbation timescale, or more precisely by the ratio of the pulse frequency,
$\omega_{j}$, and the resonant frequency, $\Omega_{\rm res}$. The right panel
of Fig. 1 shows the breathing-to-bending ratio, $f_{1,200}/f_{1,100}$, as a
function of $\omega_{\rm e}$ and $\omega_{\rm o}$, with blue (yellow) shades
indicating low (high) values. In general, the breathing-to-bending ratio rises
steeply and falls gradually with $\omega_{\rm e}$ at fixed $\omega_{\rm o}$
while the trend is reversed as a function of $\omega_{\rm o}$ at fixed
$\omega_{\rm e}$, resulting in a saddle point at $(\omega_{\rm e},\omega_{\rm
o})\approx(9,7)$. This owes to the super-exponential suppression in the
adiabatic ($\omega_{j}\ll\Omega_{\rm res}$) limit and the power-law
suppression in the impulsive ($\omega_{j}\gg\Omega_{\rm res}$) limit. Along
the $\omega_{\rm o}=\omega_{\rm e}$ line, the bending (breathing) modes
dominate in the adiabatic (impulsive) limit, as evident from the left panel of
Fig. 1. All this suggests that bending modes dominate over breathing modes
when (i) the anti-symmetric perturbation is more impulsive, i.e., evolves
faster than the symmetric one, or (ii) both symmetric and anti-symmetric
perturbations occur over comparable timescales but slower than the stellar
vertical oscillation period.

##### Persistent spirals and bars

Next we consider perturbations caused by a persistent spiral arm or bar that
grows exponentially until it saturates at a constant strength. The
corresponding temporal modulation ${\cal M}_{j}(t)$ is given by the second of
equations (22). In this case, as shown by equation (19) of Banik & van den
Bosch (2021a),

$\displaystyle{\cal P}_{nlm}^{(j)}({\bf
I},t)=\frac{\exp{\left[\gamma_{j}t\right]}\exp{\left[-im\Omega_{\rm
P}t\right]}}{\gamma_{j}+i\Omega_{\rm
res}}\left[1-\Theta(t)\right]+i\left[\frac{\gamma_{j}\exp{\left[-i(n\Omega_{z}+l\kappa+m\Omega_{\phi})t\right]}}{\Omega_{\rm
res}(\gamma_{j}+i\Omega_{\rm res})}-\frac{\exp{\left[-im\Omega_{\rm
P}t\right]}}{\Omega_{\rm res}}\right]\Theta(t).$ (31)

Up to $t=0$ when the perturber amplitude stops growing, the response from all
modes oscillates with the pattern speed $\Omega_{\rm P}$ and grows hand in
hand with the perturber. Subsequently, as the perturbation attains a steady
strength, the disk response undergoes temporary phase mixing due to the
oscillations of stars at different frequencies, giving rise to phase-spirals.
These transients, however, are quickly taken over by long term oscillations
driven at the forcing frequency $\Omega_{\rm P}$.

For a slowly growing spiral/bar, i.e., in the ‘adiabatic growth’ limit
($\gamma\to 0$), the entire disk oscillates at the driving frequency,
$\Omega_{\rm P}$, i.e.,

$\displaystyle{\cal P}_{nlm}^{(j)}({\bf I},t)$
$\displaystyle\xrightarrow{\gamma_{j}\to 0}\exp{\left[-im\Omega_{\rm
P}t\right]}\left[\pi\delta(\Omega_{\rm res})-\frac{i}{\Omega_{\rm
res}}\right].$ (32)

This has two major implications. First of all, since all stars, both resonant
and non-resonant, are driven at the pattern speed of the perturbing
spiral/bar, transient phase mixing does not occur and thus no phase-spiral
arises. Secondly, the response is dominated by the resonances, $\Omega_{\rm
res}=0$. In fact the resonant response diverges, reflecting the failure of
(standard) linear perturbation theory near resonances. The adiabatic
invariance of actions is partially broken near these resonances, causing the
stars to get trapped in librating near-resonant orbits. A proper treatment of
the near-resonant response can be performed by working with ‘slow’ and ‘fast’
action-angle variables (Tremaine & Weinberg, 1984; Lichtenberg & Lieberman,
1992; Chiba & Schönrich, 2022; Banik & van den Bosch, 2022; Hamilton et al.,
2022), which are uniquely defined for each resonance as linear combinations of
the original action-angle variables. The fast actions remain nearly invariant
while the fast angles oscillate with periods comparable to the unperturbed
orbital periods of stars. The slow action-angle variables, on the other hand,
undergo large amplitude oscillations about their resonance values over a
libration timescale that is typically much longer than the orbital periods.
For example, at co-rotation resonance ($n=l=0$), angular momentum behaves as
the slow action while the radial and vertical actions behave as the fast ones.

Based on the above discussion, we infer that phase-spirals can only be excited
in the galactic disk by transient spiral/bar perturbations whose amplitude
changes over a timescale comparable to the vertical oscillation periods of
stars. Persistent spirals or bars rotating with a fixed pattern speed cannot
give rise to phase-spirals. Rather they trigger stellar oscillations at the
pattern speed itself, which manifests in phase-space as a steadily rotating
dipole or quadrupole depending on whether the $n=1$ or $2$ mode dominates the
response. Thus, a phase-spiral is necessarily always triggered by a transient
perturbation.

####  3.1.2 Impact of collisions on the disk response

Figure 2: Timescale at which the disk response damps away due to collisional
diffusion, i.e., small-scale scatterings of stars with structures like GMCs,
is plotted as a function of $I_{z}$ ($R_{c}$) for three different values of
$R_{c}$ ($I_{z}$) as indicated, in the left (right) panel. Typically,
collisional diffusion occurs faster for smaller $I_{z}$ and smaller $R_{c}$.

In the above section we discussed the characteristics of the disk response in
the absence of collisions. However, in a real galaxy like the MW disk, small-
scale collisionality can potentially damp away any coherent response to a
perturbation. Collisional diffusion arises not from star-star collisions,
which is typically negligible, but from gravitational scattering with other
objects, such as GMCs, DM substructure, etc. As discussed in Section 2.2, the
impact of collisional diffusion is mainly captured by the diffusion
coefficients $D_{z}$ and $D_{R}$. Following Tremaine et al. (2022) we assume
that the disk stars have gained their mean vertical and radial actions over
the age of the disk, $T_{\rm disk}=10\>{\rm Gyr}$, due to collisional heating,
which implies that $D_{a}=\left<I_{a}\right>/T_{\rm disk}$ where $a$ is either
$z$ or $R$ and $\left<I_{a}\right>=\int{\rm d}I_{a}\,I_{a}\,f_{0}\,/\int{\rm
d}I_{a}\,f_{0}$.

Figure 3: MW disk response to transient bars/2-armed spirals with Gaussian
temporal modulation of pulse frequency, $\omega_{\rm o}=\omega_{\rm
e}=0.5\,\sigma_{z,\odot}/h_{z}$: the amplitude of the disk response,
$f_{1n00}/f_{0}$, is plotted as a function of time. The rows and columns
respectively denote different values of $R_{c}$ and $I_{z}$ as indicated. Blue
and red lines indicate the $n=1$ and $2$ modes, while the solid and dashed
lines respectively denote the cases with and without collisional diffusion
(due to interactions of stars with structures like GMCs). The disk response
initially rises and falls hand in hand with the perturbing pulse (indicated by
the grey dotted line), before saturating to a steady state in the
collisionless case and undergoing super-exponential damping in the collisional
case. Note that collisional damping is faster for smaller $I_{z}$, smaller
$R_{c}$ and larger $n$ modes.

For a transient bar/spiral with pulse frequency $\omega_{j}$, ${\cal
P}^{(j)}_{nlm}$ is given by equation (25). In the impulsive limit
($\omega_{j}\to\infty$), we have that ${\cal
M}_{j}(t-\tau)\to\omega_{j}\delta(t-\tau)$. Upon absorbing $\omega_{j}$ in the
prefactor, the expression for ${\cal P}^{(j)}_{nlm}$ then simplifies to (c.f.
Appendix A of Tremaine et al., 2022)

$\displaystyle{\cal P}^{(j)}_{nlm}({\bf I},t)$
$\displaystyle\approx\Theta(t)\,\exp{\left[-i\,(n\Omega_{z}+l\kappa+m\Omega_{\phi})\,t\right]}\,\exp{\left[-\left(\frac{n^{2}D_{z}}{4I_{z}}+\frac{l^{2}D_{R}}{4I_{R}}\right)t\right]}\,\exp{\left[-\frac{{\left(n\Omega_{z1}\right)}^{2}D_{z}I_{z}}{3}t^{3}\right]}.$
(33)

This demonstrates that, in the impulsive limit, the disk response
instantaneously grows and spawns phase-spirals whose amplitude decays due to
collisional diffusion, manifest in the exponential damping terms. The first
and second exponential factors, respectively, characterize the diffusion in
vertical angle and action, which occur over the following timescales:

$\displaystyle\tau_{\rm D}^{({\rm w})}$
$\displaystyle={\left[\frac{n^{2}D_{z}}{4I_{z}}+\frac{l^{2}D_{R}}{4I_{R}}\right]}^{-1},\;\;\;\;\tau_{\rm
D}^{({\rm
I})}={\left[\frac{3}{{\left(n\Omega_{z1}\right)}^{2}D_{z}I_{z}}\right]}^{1/3}.$
(34)

Of these, the timescale for the diffusion in angles, $\tau_{\rm D}^{({\rm
w})}$, typically exceeds that for the diffusion in actions, $\tau_{\rm
D}^{({\rm I})}$, by at least an order of magnitude, implying that angle
diffusion is negligible. Hence collisional diffusion mainly causes the
abatement of action gradients in the phase-space structure of the response
(arising from the action dependence of the frequencies, i.e., $\Omega_{z1}\neq
0$). The left (right) panel of Fig. 2 plots the diffusion timescale,
$\tau_{\rm D}^{({\rm I})}$, as a function of $I_{z}$ ($R_{\rm c}$) for three
different values of $R_{\rm c}$ ($I_{z}$) as indicated. Note that $\tau_{\rm
D}^{({\rm I})}$ diverges in the small $I_{z}$ limit, attains a minimum around
$I_{z}\sim 0.2-0.5\,h_{z}\sigma_{z,\odot}$, and increases as $I_{z}^{\beta}$
with $\beta<1$ at large $I_{z}$. As a function of $R_{\rm c}$, $\tau_{\rm
D}^{({\rm I})}$ shows an approximately exponential rise. This owes to the fact
that $\left<D_{z}\right>\sim h_{z}\,\sigma_{z}(R_{\rm c})/T_{\rm
disk}\sim\exp{\left[-R_{\rm c}/2h_{R}\right]}/T_{\rm disk}$ for the 3MN
profile adopted for the MW disk. At $R_{\rm c}=\>{R_{\odot}}=8\>{\rm kpc}$,
$\tau_{\rm D}^{(I)}\sim 0.6-0.7\>{\rm Gyr}$, in agreement with Tremaine et al.
(2022). Hence, we see that collisional diffusion in action space is fairly
efficient, and thus that phase-spirals are short-lived features.

Fig. 3 plots the amplitude of the disk response (for $I_{R}=0$) to a transient
spiral of pulse frequency, $\omega_{\rm o}=\omega_{\rm
e}=0.5\sigma_{z,\odot}/h_{z}$, computed using equations (20), (24) and (25),
as a function of time. Dashed and solid lines show the results with and
without collisional diffusion, respectively. The rows correspond to different
values of $R_{\rm c}$ while the columns denote different values of
$I_{z}/(h_{z}\sigma_{z,\odot})$ as indicated. The blue and red lines denote
the response for the $(n,l,m)=(1,0,0)$ and $(2,0,0)$ modes, respectively, and
the dotted grey line represents the Gaussian pulse strength. The response for
both bending and breathing modes initially grows hand in hand with the
perturbing pulse. Following the point of maximum pulse strength, the response
follows the decaying pulse strength before saturating to the steady state
amplitude given in equation (27) in the collisionless limit. In the presence
of collisional diffusion, however, the response continues to damp out as
$\sim\exp{[-(t/\tau_{\rm D}^{({\rm I})})^{3}]}$ after temporarily saturating
at the collisionless steady state. Note that the collisional damping is faster
for smaller $R_{\rm c}$ and smaller $I_{z}$. In addition, $n=2$ breathing
modes damp out faster than the $n=1$ bending modes due to the $n^{-2/3}$
dependence of $\tau_{\rm D}^{({\rm I})}$.

To summarize, we have shown that phase-spirals can be triggered by impulsive
perturbations resulting from transient spiral arms or bars, but are subject to
super-exponential damping due to collisional diffusion that is likely to be
dominated by scattering against GMCs. This collisional damping is more
efficient in the inner disk, for stars with smaller $I_{z}$, and for modes of
larger $n$.

##  4 Disk response to satellite encounter

In addition to the spiral arm/bar perturbations considered above, we also
consider disk perturbations triggered by encounters with a satellite galaxy.
For the sake of brevity, we only compute the disk response in the
collisionless limit. In the case of impulsive encounters, the impact of
collisional diffusion is simply expressed by multiplying the collisionless
response by the collisional damping factor $\exp[-(t/\tau_{\rm D}^{({\rm
I})})^{3}]$, with $\tau_{\rm D}^{({\rm I})}$ given by equation (34).

For simplicity, we assume that the satellite is moving with uniform velocity
$v_{\rm P}$ along a straight line, impacting the disk at a galactocentric
distance $R_{\rm d}$ with an arbitrary orientation, specified by the angles,
$\theta_{\rm P}$ and $\phi_{\rm P}$, which are respectively defined as the
angles between ${{\bf v}}_{\rm{{\bf P}}}$ and the $z$-axis, and between the
projection of ${{\bf v}}_{\rm{{\bf P}}}$ on the mid-plane and the $x$-axis
(see Fig. 4). Thus the position vector of the satellite with respect to the
galactic center can be written as

$\displaystyle{{\bf r}}_{\rm{{\bf P}}}=(R_{\rm d}+v_{\rm P}\sin{\theta_{\rm
P}}\cos{\phi_{\rm P}}\,t)\,\hat{{\bf x}}+v_{\rm P}\sin{\theta_{\rm
P}}\sin{\phi_{\rm P}}\,t\,\hat{{\bf y}}+v_{\rm P}\cos{\theta_{\rm
P}}\,t\,\hat{{\bf z}},$ (35)

while that of a star is given by

$\displaystyle{\bf r}=R(\cos{\phi}\,\hat{{\bf x}}+\sin{\phi}\,\hat{{\bf
y}})+z\,\hat{{\bf z}}.$ (36)

We consider the satellite to be a Plummer sphere of mass $M_{\rm P}$ and size
$\varepsilon$, such that its gravitational potential at location ${\bf r}$ is
given by

$\displaystyle\Phi_{\rm P}$ $\displaystyle=GM_{\rm
P}\left[-\frac{1}{\sqrt{{\left|{\bf r}-{{\bf r}}_{\rm{{\bf
P}}}\right|}^{2}+\varepsilon^{2}}}+\frac{{\bf r}\cdot{{\bf r}}_{\rm{{\bf
P}}}}{{\left({{\bf r}}_{\rm{{\bf
P}}}^{2}+\varepsilon^{2}\right)}^{3/2}}\right]\,.$ (37)

Here the first term is the ‘direct’ term and the second is the ‘indirect’ term
that accounts for the reflex motion of the disk and the fact that the disk
center is accelerated by the satellite and is thus non-inertial. Typically,
the first one dominates over the second.

In order to compute the disk response to this external perturbation, we need
to compute its Fourier coefficients, which is challenging. Rather, we first
evaluate the $\tau$-integral in Equation (20), setting $t_{\rm i}\to-\infty$,
and then compute the Fourier transform of the result, as worked out in
Appendix C.1. For $I_{R}\approx 0$ (this is justified since we adopt the
radial epicyclic approximation in this paper), this yields a modal response,
$f_{1nlm}$ (Equation [20]), with ${\cal I}_{nlm}({\bf I},t)$ given by

$\displaystyle{\cal I}_{nlm}({\bf I},t)\approx-\frac{2GM_{\rm P}}{v_{\rm
P}}\exp{\left[-i\Omega
t\right]}\times\exp{\left[-i\frac{\Omega\sin{\theta_{\rm P}}\cos{\phi_{\rm
P}}}{v_{\rm P}}R_{\rm d}\right]}$
$\displaystyle\times\frac{1}{{\left(2\pi\right)}^{2}}\int_{0}^{2\pi}dw_{z}\exp{\left[-inw_{z}\right]}\exp{\left[i\frac{\Omega\cos{\theta_{\rm
P}}}{v_{\rm
P}}z\right]}\int_{0}^{2\pi}d\phi\,\exp{\left[-im\phi\right]}\,\exp{\left[i\frac{\Omega\sin{\theta_{\rm
P}}\cos{\left(\phi-\phi_{\rm P}\right)}}{v_{\rm P}}R_{\rm c}\right]}$
$\displaystyle\times K_{0i}\left(\frac{\Omega\sqrt{{\cal R}^{2}_{\rm
c}+\varepsilon^{2}}}{v_{\rm P}},\frac{v_{\rm P}t-{\cal S}_{\rm c}}{\sqrt{{\cal
R}^{2}_{\rm c}+\varepsilon^{2}}}\right),$ (38)

where $\Omega$ is given by

$\displaystyle\Omega=n\Omega_{z}+l\kappa+m\Omega_{\phi}.$ (39)

Here ${\cal R}_{\rm c}={\cal R}(R_{\rm c})$ and ${\cal S}_{\rm c}={\cal
S}(R_{\rm c})$ with ${\cal R}$ and ${\cal S}$ given by equation (C5). $K_{0i}$
is given by equation (C3), which asymptotes to the modified Bessel function of
the second kind, $K_{0}\left(\left|\Omega\right|\sqrt{{\cal R}^{2}_{\rm
c}+\varepsilon^{2}}/v_{\rm P}\right)$, in the large time limit. A more precise
expression for ${\cal I}_{nlm}$ that is valid for higher values of $I_{R}$ is
given by equation (C7) of Appendix C.1.

Figure 4: Illustration of the geometry of a satellite galaxy with mass $M_{\rm
P}$ impacting a disk galaxy with uniform velocity $v_{\rm P}$ along a straight
line. The impact occurs at a galactocentric distance $R_{\rm d}$. The
orientation of ${{\bf v}}_{\rm{{\bf P}}}$ is specified by $\theta_{\rm P}$,
the angle between ${{\bf v}}_{\rm{{\bf P}}}$ and the $z$-axis, and $\phi_{\rm
P}$, the angle between the projection of ${{\bf v}}_{\rm{{\bf P}}}$ on the
mid-plane and the $x$-axis.

The expression for ${\cal I}_{nlm}$ given in equation (38) exhibits several
key features of the disk response to satellite encounters. The
$\exp{\left[-i\Omega t\right]}$ factor encodes the phase mixing of the
response due to oscillations at different frequencies, giving rise to phase-
spirals. The $\exp{\left[i\left(\Omega\cos{\theta_{\rm P}}/v_{\rm
P}\right)z\right]}$ and $\exp{\left[i\left(\Omega\sin{\theta_{\rm
P}}\cos{(\phi-\phi_{\rm P})}/v_{\rm P}\right)R_{\rm c}\right]}$ factors
respectively indicate that the satellite induces wave-like perturbations in
the disk with two characteristic wave-numbers: the vertical wave-number,
$k_{z}\approx\Omega\cos{\theta_{\rm P}}/v_{\rm P}$ and the horizontal wave-
number, $k_{R}\approx\Omega\sin{\theta_{\rm P}}/v_{\rm P}$. Therefore, the
disk response will be vertically (horizontally) stratified in case of a
perpendicular (planar) impact of the satellite. As shown in Appendix C.2,
expressions (20) and (38), which are complicated to compute, yield the correct
response in the impulsive limit of a satellite having a face-on, perpendicular
encounter through the center of the disk.

###  4.1 Asymptotic behaviour of the response

It is instructive to study the two extreme cases of encounter speed, the
impulsive limit (large $v_{\rm P}$) and the adiabatic limit (small $v_{\rm
P}$). Using the asymptotic form of the $K_{0}$ Bessel function that appears in
equation (38), we obtain the following approximate asymptotic behaviour of
$f_{1nlm}$ at large time:

$\displaystyle f_{1nlm}$ $\displaystyle\sim\frac{GM_{\rm P}}{v_{\rm
P}}\exp{\left[-i\Omega t\right]}\times\begin{cases}1,&v_{\rm P}\to\infty\\\
\\\ \sqrt{v_{\rm P}/\Omega b}\,\exp{\left[-\Omega b/v_{\rm P}\right]},&v_{\rm
P}\to 0,\end{cases}$ (40)

where $b$ is the impact parameter of the encounter, defined as the
perpendicular distance of the nearest star on the mid-plane from the
satellite’s (straight) orbit, and expressed as

$\displaystyle b=\left|R_{\rm d}-R_{\rm
c}\right|\,\sqrt{1-\sin^{2}{\theta_{\rm P}}\cos^{2}{\phi_{\rm P}}}.$ (41)

It is clear from these limits that the disk response is most pronounced for
intermediate velocities, $v_{\rm P}\sim\Omega b$. For impulsive encounters,
the response is suppressed as a power law in $v_{\rm P}$, whereas in the
adiabatic limit the response is exponentially suppressed, except at
resonances, $\Omega=n\Omega_{z}+l\kappa+m\Omega_{\phi}=0$. In this limit, far
from the resonances, the perturbation timescale, $b/v_{\rm P}$, is much larger
than $\Omega^{-1}$, and the net response is washed away due to many
oscillations during the perturbation (i.e., the actions are adiabatically
invariant), a phenomenon known as adiabatic shielding (Weinberg, 1994a, b;
Gnedin & Ostriker, 1999).

Figure 5: Steady state MW disk response to satellite encounter in the
collisionless limit: bending mode strength, $f_{1,n=1}/f_{0}$ (upper panel),
and the corresponding breathing vs bending ratio, $f_{1,n=2}/f_{1,n=1}$ (lower
panel) for the $(l,m)=(0,0)$ modes, in the Solar neighborhood for the MW
satellites, as a function of the disk crossing time, $t_{\rm cross}$, in Gyr,
where $t_{\rm cross}=0$ marks today. The previous two and the next impacts are
shown. Here we consider $I_{z}=h_{z}\sigma_{z,\odot}$, with fiducial MW
parameters, and marginalize over $I_{R}$. The effect of the (non-responsive)
ambient DM halo on the stellar frequencies is taken into account. The
estimates of $t_{\rm cross}$ are very sensitive to the detailed potential of
the MW system, while the response estimates are fairly robust (see text for
details). In the upper panel, the region with bending mode response,
$f_{1,n=1}/f_{0}<10^{-4}$, has been grey-scaled, indicating that the response
from the satellites in this region is far too weak and adiabatic to be
detected by Gaia. Note that the response is dominated by that due to Sgr,
followed by Hercules, Leo II, Segue 2 and the Large Magellanic Cloud (LMC).
Also note that the previous two and next impacts of all the satellites excite
bending modes in the Solar neighborhood.

Figure 6: Steady state MW disk response to satellite encounter in the
collisionless limit: each panel shows the behaviour of the disk response
amplitude, $f_{1,n00}/f_{0}$ (evaluated using equations [20] and [C7]) and
marginalized over $I_{R}$), as a function of the impact velocity, $v_{\rm P}$,
in the Solar neighborhood, i.e., $R_{\rm c}=\>{R_{\odot}}=8\>{\rm kpc}$, in
presence of an ambient DM halo. The left and right columns respectively
indicate the response for the $n=1$ bending and $n=2$ breathing modes. The
top, middle and bottom rows show the same for different values of $I_{z}$ (in
units of $I_{z,\odot}$), $\theta_{\rm P}$ and $\phi_{\rm P}$ respectively as
indicated, with the fiducial parameters corresponding to $I_{z,\odot}$ and the
parameters for Sgr impact, the response amplitude for which is indicated by
the red circle. Note that the response is suppressed as $v_{\rm P}^{-1}$ in
the impulsive (large $v_{\rm P}$) limit but exponentially suppressed in the
adiabatic (small $v_{\rm P}$) regime, and peaks at an intermediate velocity,
$v_{\rm P}\sim 2-3\,v_{\rm circ}(\>{R_{\odot}})$ (which is very similar to the
encounter speed of Sgr). The peak of the response shifts to smaller $v_{\rm
P}$ for larger $I_{z}$, since $\Omega_{z}$ decreases with $I_{z}$. The
response depends only very weakly on $\phi_{\rm P}$ but is quite sensitive to
$\theta_{\rm P}$; more planar encounters, i.e., increasing $\theta_{\rm P}$
triggers stronger responses.

###  4.2 Response of the MW disk to satellites

The MW halo harbors several fairly massive satellite galaxies that repeatedly
perturb the MW disk. Here we use existing data on the phase-space coordinates
of those MW satellites to compute the disk response of satellite encounters
that occurred in the past few hundred Myr, which are those for which we may
expect phase-spirals that were triggered to have survived to the present day.

To compute the disk response to the MW satellites, we proceed as follows. As
in Paper I, we adopt the galactocentric coordinates and velocities computed
and documented by Riley et al. (2019) (table A.2, see also Li et al., 2020)
and Vasiliev & Belokurov (2020) as initial conditions for the MW satellites.
We then simulate their orbits in the combined gravitational potential of the
MW halo, disk plus bulge222The bulge is modelled as a spherical Hernquist
(1990) profile with mass $M_{b}=6.5\times 10^{9}\>{\rm M_{\odot}}$ and scale
radius $r_{b}=0.6$ kpc. using a second order leap-frog integrator. For each
individual orbit, we record the times, $t_{\rm cross}$, and the galactocentric
radii, $R_{\rm d}$, corresponding to disk crossings. We also register the
corresponding impact velocities, $v_{\rm
P}=\sqrt{v^{2}_{z}+v^{2}_{R}+v^{2}_{\phi}}$, and the angles of impact,
$\theta_{\rm P}=\cos^{-1}{(v_{z}/v_{\rm P})}$ and $\phi_{\rm
P}=\tan^{-1}{(v_{\phi}/v_{R})}$. We substitute these quantities in equation
(C7) and compute the disk response (integrated over $I_{R}$) following the
satellite encounter, using equations (20) and (38). Results are summarized in
Table 1 of Appendix C.1. Fig. 5 plots the amplitude of the Solar neighborhood
(for which $R_{\rm c}(L_{z})=\>{R_{\odot}}=8\>{\rm kpc}$) bending mode
response, $f_{1,n=1}/f_{0}$ (top panel), and breathing-to-bending ratio,
$f_{1,n=2}/f_{1,n=1}$ (bottom panel), as a function of $t_{\rm cross}$. Here
we only show the responses for $(l,m)=(0,0)$ modes, and consider stars with
$I_{z}=I_{z,\odot}=9.2\>{\rm kpc}\>{\rm km}\,{\rm s}^{-1}$.

It is noteworthy that the responses in the realistic MW disk computed here are
$\sim 1-2$ orders of magnitude larger than those evaluated for the isothermal
slab model shown in Fig. 7 of Paper I. This owes to the reduced damping of the
phase-spiral amplitude due to lateral mixing, which is more pronounced in the
isothermal slab with unconstrained lateral velocities than in the realistic
disk with constrained, ordered motion. From the lower panel of Fig. 5 it is
evident that, as in the isothermal slab case, almost all satellites trigger a
bending mode response in the Solar neighborhood, resulting in a one-armed
phase-spiral in qualitative agreement with the Gaia snail. However, as is
evident from the upper panel, only five of the satellites trigger a detectable
response in the disk, with $f_{1,n=1}/f_{0}>\delta_{\rm min}\equiv 10^{-4}$
(see Appendix C of Paper I for a derivation of this approximate detectability
criterion for Gaia). The response to encounters with the other satellites is
weak either because they have too low mass or because the encounter with
respect to the Sun is too slow and adiabatically suppressed. Sgr excites the
strongest response by far; its bending mode response, $f_{1,n=1}/f_{0}$, is at
least $1-2$ orders of magnitude above that for any other satellite. Its
penultimate disk crossing, about the same time as its last pericentric passage
$\sim 1\>{\rm Gyr}$ ago, triggered a strong response of $f_{1,n=1}/f_{0}\sim
0.3$ in the Solar neighborhood. For comparison, the response from its last
disk crossing, which nearly coincides with its last apocentric passage about
$350\>{\rm Myr}$ ago, triggered a very weak, adiabatically suppressed response
($\sim 5\times 10^{-8}$) that falls below the lower limit of Fig. 5. Its next
disk crossing in about $30\>{\rm Myr}$ is estimated to trigger a strong
response with $f_{1,n=1}/f_{0}\sim 0.1$. Besides Sgr, the satellites that
excite a detectable response, $f_{1,n=1}/f_{0}>\delta_{\rm min}$ are Hercules,
Segue 2, Leo II and the LMC. The imminent crossing of LMC is estimated to
trigger $f_{1,n=1}/f_{0}\sim 2\times 10^{-2}$, which is an order of magnitude
below Sgr. Only for $I_{z}/I_{z,\odot}\gtrsim 4.5$ ($z_{\rm max}\gtrsim
3.4h_{z}$), the LMC response dominates over Sgr. This exercise therefore
suggests that Sgr is the leading contender, among the MW satellites considered
here, for triggering the Gaia snail in the Solar neighborhood, in agreement
with several previous studies (Antoja et al., 2018; Binney & Schönrich, 2018;
Laporte et al., 2018, 2019; Darling & Widrow, 2019a; Bland-Hawthorn et al.,
2019; Hunt et al., 2021; Bland-Hawthorn & Tepper-García, 2021; Bennett et al.,
2022).

We caution that the estimates of the disk response computed above ignore
dynamical friction. Moreover, the disk crossing times are sensitive to the
satellite orbits and therefore to the detailed MW potential and the current
phase-space coordinates of the satellites. For example, a heavier MW model
with a total mass of $1.5\times 10^{12}\>{\rm M_{\odot}}$ leaves the relative
amplitudes of the satellite responses (in the collisionless limit) nearly
unchanged, but makes the satellites more bound, bringing most of the disk
crossing times closer to the present day. In particular, the last disk passage
of Sgr that triggers a significant response now occurs $\sim 600\>{\rm Myr}$
ago (as opposed to $1\>{\rm Gyr}$ ago in the fiducial case) which is closer to
the winding time of $\sim 500\>{\rm Myr}$ inferred from the phase-spiral
observed in the Solar neighborhood (Bland-Hawthorn et al., 2019).

In this section, we have computed the responses in the collisionless limit. In
reality, collisional diffusion due to interactions of stars with GMCs, etc.
would damp away the response super-exponentially over a timescale that is
$\sim 0.6-0.7\>{\rm Gyr}$ in the Solar neighborhood (see section 3.1.2). This
would almost completely wash away the response to any satellite encounter that
occurred $\gtrsim 1\>{\rm Gyr}$ ago. For example, the present day response to
the last pericentric passage of Leo II that occurred $\sim 1.8\,\>{\rm Gyr}$
ago would be completely erased. If the last disk crossing of Sgr that induced
a strong response occurred $\sim 1\>{\rm Gyr}$ ago as in the fiducial MW
model, the response would have been damped out by $\sim 2$ orders of magnitude
by today, deeming Sgr unlikely to be the agent behind the Gaia snail. However,
as discussed above, the disk crossing times are sensitive to the satellite
orbits. The heavier MW model with a total mass of $1.5\times 10^{12}\>{\rm
M_{\odot}}$ implies a Sgr crossing time of $\sim 0.6\>{\rm Gyr}$ instead of
$1\>{\rm Gyr}$. In this case the response would only have been damped by a
factor of $\sim 0.4$. Therefore, the collisionality argument suggests that if
the Gaia snail was indeed triggered by Sgr, the impact causing it must have
happened within $\sim 0.6-0.7\>{\rm Gyr}$ from the present day.

###  4.3 Exploring parameter space

Figure 7: Steady state MW disk response to satellite encounter in the
collisionless limit: each panel shows the behaviour of the disk response
amplitude, $f_{1,nlm}/f_{0}$ (marginalized over $I_{R}$), as a function of the
impact velocity, $v_{\rm P}$, in the Solar neighborhood, in presence of an
ambient DM halo. Different lines correspond to different $m$ modes as
indicated. The top and bottom rows show the response for $l=0$ and $1$ while
the left and right columns indicate it for the $n=1$ bending and $n=2$
breathing modes. The fiducial parameters correspond to $I_{z}=I_{z,\odot}$ and
the parameters for Sgr impact, the response amplitudes for which are indicated
by the red circles in each panel. The response is dominated by the
$(n,l,m)=(1,0,-2)$ mode or the two-armed warp at small $v_{\rm P}$ and the
$(2,0,-2)$ mode or the two-armed spiral at large $v_{\rm P}$. Typically, the
$m=-2$ and $-1$ responses dominate over $m=0,1$ and $2$, while the $l=0$
response is more pronounced than $l=1$.

Figure 8: MW disk response to satellite encounter: breathing-to-bending ratio
or the relative strength of the $n=2$ and $n=1$ modes of disk response to a
Sgr-like impact is plotted as a function of the impact velocity, $v_{\rm P}$,
at $R_{\rm c}=\>{R_{\odot}}=8\>{\rm kpc}$ and $R_{\rm
c}=1.5\>{R_{\odot}}=12\>{\rm kpc}$ shown in the left and right columns
respectively, for the $(l,m)=(0,-2)$ mode which typically dominates the
response. Different lines correspond to different values of $\theta_{\rm P}$
as indicated. We consider $I_{z}=I_{z,\odot}$ and the fiducial parameters to
correspond to those for Sgr encounter, for which the breathing-to-bending
ratio is denoted by the red circle. Bending modes dominate over breathing
modes at small $v_{\rm P}$ and vice versa at large $v_{\rm P}$. Breathing
modes are relatively more pronounced than bending modes in the outer disk,
closer to the Sgr impact radius, $R_{\rm d}=17\>{\rm kpc}$. More planar
(perpendicular) encounters trigger larger breathing-to-bending ratios farther
away from (closer to) the impact radius.

Having computed the MW disk response to its satellites, we now investigate the
sensitivity of the response to the various encounter parameters. In Fig. 6 we
plot the amplitude of the Solar neighborhood response, $f_{1,nlm}/f_{0}$
(marginalized over $I_{R}$), as a function of the impact velocity, $v_{\rm P}$
(in units of the circular velocity at $R_{\rm c}=\>{R_{\odot}}$), for the
$(n,l,m)=(1,0,0)$ bending and $(n,l,m)=(2,0,0)$ breathing modes, shown in the
left and right columns respectively. The top, middle and bottom rows show the
results for varying $I_{z}$, $\theta_{\rm P}$ and $\phi_{\rm P}$ respectively,
assuming the fiducial parameters to be those for Sgr (mass $M_{\rm
P}=10^{9}\>{\rm M_{\odot}}$, scale radius $\varepsilon=1.6$ kpc) during its
penultimate disk crossing (most relevant for the Gaia snail), i.e., impact
radius $R_{\rm d}=17\>{\rm kpc}$, impact velocity $v_{\rm P}=340$ km/s, and
angles of impact, $\theta_{\rm P}={21}^{\circ}$ and $\phi_{\rm
P}={150}^{\circ}$. In Fig. 7 we plot the bending and breathing mode response
amplitudes (in the Solar neighborhood) as a function of $v_{\rm P}$ for
different $(l,m)$ modes, with the fiducial parameters again corresponding to
Sgr. The left and right columns respectively indicate the $n=1$ bending and
$n=2$ breathing modes, while the top and bottom rows correspond to $l=1$ and
$l=2$ respectively. The different lines in each panel denote the responses for
$m=-2,-1,0,1$ and $2$. Fig. 8 shows the ratio of the bending and breathing
response amplitudes as a function of $v_{\rm P}$ for the dominant mode
$(l,m)=(0,-2)$. Different lines indicate breathing-to-bending ratios for
different values of $\theta_{\rm P}$, while the left and right columns
respectively indicate the ratios observed at $R_{\rm c}=8$ and $12\>{\rm
kpc}$.

From Figs. 6 and 7 it is evident that, as shown in equation (40), the disk
response is suppressed like a power law ($\sim v^{-1}_{\rm P}$) in the high
velocity/impulsive limit and exponentially ($\sim\exp{\left[-\Omega b/v_{\rm
P}\right]}$) suppressed in the low velocity/adiabatic limit. The response is
the strongest for intermediate velocities, $v_{\rm P}\sim 2-3\,v_{\rm
circ}(\>{R_{\odot}})$, where the time periods of the vertical, radial and
azimuthal oscillations of the stars are nearly commensurate with the encounter
timescale, $\sqrt{b^{2}+\varepsilon^{2}}/v_{\rm P}$. The $v^{-1}_{\rm P}$ and
$K_{0i}$ factors in equation (38) conspire to provide the near-resonance
condition for maximum response,

$\displaystyle n\Omega_{z}+l\kappa+m\Omega_{\phi}\approx\frac{0.6\,v_{\rm
P}}{\sqrt{b^{2}+\varepsilon^{2}}},$ (42)

where $b$ is the impact parameter of the encounter, given by equation (41).
From the top panels of Fig. 6, it is clear that the peak response shifts to
smaller $v_{\rm P}$ with increasing $I_{z}$. This is easy to understand from
the fact that the corresponding vertical frequency, $\Omega_{z}$, decreases
with increasing $I_{z}$, making the encounter more impulsive for larger
actions. The middle and bottom panels show that the response depends strongly
on the polar angle of the encounter, $\theta_{\rm P}$, but very mildly on the
azimuthal angle, $\phi_{\rm P}$. Moreover, the middle panels indicate that
more planar encounters (larger $\theta_{\rm P}$) induce stronger responses.

The in-plane structure of the disk response depends on the relative
contribution of the different $(l,m)$ modes. From Fig. 7 it is evident that a
typical Sgr-like encounter predominantly excites $(l,m)=(0,-1)$ and
$(l,m)=(0,-2)$ in the Solar neighborhood. The dominant mode for slower
encounters is $(n,l,m)=(1,0,-2)$ while that for faster ones is
$(n,l,m)=(2,0,-2)$. Since $f_{1,nlm}/f_{0}\gtrsim 1$ in these cases, the
response to the impact by Sgr is in fact non-linear in the Solar neighborhood.
Either way, a satellite encounter is typically found to excite strong $m=-2$
modes, i.e., $2$-armed warps ($n=1$) and spirals ($n=2$). This is due to a
quadrupolar tidal distortion of the disk by the satellite, which manifests as
a stretching of the disk in the direction of the impact and a compression
perpendicular to it.

Fig. 8 elucidates that the bending mode response dominates for slower
encounters, i.e., smaller $v_{\rm P}$, and at guiding radii far from the
impact radius, $R_{\rm d}$. More planar impacts trigger larger breathing-to-
bending ratios farther away from the impact radius while this trend reverses
closer to it. This is because more planar encounters cause more vertically
symmetric perturbations farther away from the impact radius. The predominance
of bending modes for low $v_{\rm P}$ encounters while that of breathing modes
for high $v_{\rm P}$ ones has been observed by Widrow et al. (2014) and Hunt
et al. (2021) in their N-body simulations of satellite-disk encounters. As
demonstrated by Widrow et al. (2014), slower encounters provide energy to the
stars near one of the vertical turning points while drain energy from those
near the other turning point, thereby driving bending wave perturbations that
are asymmetric about the mid-plane. On the other hand, fast satellite passages
are impulsive and impart energy to the stars near both the turning points,
thus triggering symmetric breathing waves.

The predominance of breathing (bending) modes closer to (farther away from)
the impact radius is qualitatively similar to the observation by Hunt et al.
(2021) in their simulations of MW-Sgr encounter that the outer part of the MW
disk which is closer to the impact radius shows a preponderance of two-armed
phase-spirals or breathing modes. This can be understood within the framework
of our formalism by noting that the impact parameter, $b$, and therefore the
encounter timescale $\sim\sqrt{b^{2}+\varepsilon^{2}}/v_{\rm P}$ decreases
with increasing proximity to the point of impact; hence the impact is faster
than the vertical oscillations of stars near the point of impact, driving
stronger breathing mode perturbations. However, contrary to these predictions
for the MW-Sgr encounter, Hunt et al. (2022), using Gaia DR3 data, revealed
two-armed phase-spirals, and therefore breathing modes, in the inner disk
($R_{\rm c}\sim 6-7\>{\rm kpc}$). Our analysis suggests that none of the MW
satellites could have caused this. Using N-body simulations of an isolated MW
system, Hunt et al. (2022) suggested that a transient spiral arm or bar could
be a potential trigger for breathing modes in the inner disk. However, such a
transient perturbation would have to be sufficiently impulsive, i.e., occur
over a timescale that is comparable to or smaller than the vertical
oscillation timescale in the inner disk (see section 3.1.1), in order to
produce two-armed phase-spirals with density contrast as strong as in the
data. Such short timescales are unlikely to arise from the secular evolution
of the disk alone and may instead require forcing of the inner disk by
perturbations in the MW halo. Another possible trigger of this feature is the
recent passage of dark satellite(s) through the inner disk. The true origin of
this feature is however unclear. Hence, we conclude that the presence of two-
armed phase-spirals in the inner disk is rather unexpected, and that its
origin poses an intriguing conundrum.

##  5 phase-spirals and the Galactic potential

Thus far we mainly focused on how the nature of the perturbation dictates the
vertical (i.e., bending and breathing modes) as well as the in-plane (various
$(l,m)$ modes) structure of the disk response. However, the detailed
structure, in particular the winding, of the phase-spiral not only depends on
the triggering agent but also holds crucial information about the underlying
potential in which the stars move, and can thus be used to constrain the
potential of the combined disk plus halo system (see also Widmark et al.,
2022a, b).

The winding of the vertical phase-spiral can be characterized by the pitch-
angle, $\phi_{\rm I}$, along the ridge of maximum density. It is defined as
the angle between the azimuthal direction and the tangent to the line of
constant density (Binney & Tremaine, 1987). It is related to the local
dependence of the vertical frequency on the vertical action according to:

$\displaystyle\phi_{\rm I}=\cot^{-1}{\left[\left|I_{z}\frac{{\rm
d}\Omega_{z}}{{\rm d}I_{z}}\right|t\right]}=\cot^{-1}{\left[\left|\frac{{\rm
d}\Omega_{z}}{{\rm d}\ln{I_{z}}}\right|t\right]}.$ (43)

Following a perturbation, the pitch angle decreases with time, asymptoting
towards zero, as the spiral winds up as a consequence of the ongoing phase
mixing. Based on the above expression for $\phi_{\rm I}$, we can define the
following timescale of phase mixing:

$\displaystyle\tau_{\phi}=\left|\frac{{\rm d}\ln{I_{z}}}{{\rm
d}\Omega_{z}}\right|\,.$ (44)

This timescale, which determines the rate of winding of the spiral, is a
function of both the guiding radius, $R_{\rm c}$, and the action, $I_{z}$, and
is ultimately dictated by the (unperturbed) potential of the disk+halo system,
which sets ${\rm d}\Omega_{z}/{\rm d}I_{z}$. Hence, the detailed shape of the
phase-spiral at a given location in the disk is sensitive to the local
disk+halo potential, thereby opening up interesting avenues for constraining
the detailed potential of the MW by examining phase-spirals throughout the
disk.

Figure 9: Impact of DM halo on vertical phase mixing: the panels from left to
right respectively indicate the vertical frequency, $\Omega_{z}$ (units of
$\sigma_{z,\odot}/h_{z}$), the vertical phase mixing timescale, $\tau_{\phi}$
(given by equation [44]), and the $w_{z}=0$ cuts of the phase-spirals shown in
Fig. 10 as a function of the vertical action, $I_{z}$ (units of
$h_{z}\sigma_{z,\odot}$). The solid and dashed red lines denote the cases with
and without a halo for $R_{\rm c}=\>{R_{\odot}}=8\>{\rm kpc}$ while the dot-
dashed and dotted blue lines show the same for $R_{\rm c}=12\>{\rm kpc}$. The
vertical dashed line indicates roughly the maximum $I_{z}$ for which a phase-
spiral is discernible in the Gaia data. Note that phase mixing occurs the
fastest for $I_{z}\sim 1$ and that the inner disk phase mixes faster than the
outer disk. Also note that the presence of a DM halo increases $\Omega_{z}$ as
well as $\tau_{\phi}$, leading to slower phase mixing and therefore slower
wrapping of the phase-spiral. This effect is more pronounced in the outer
disk.

Figure 10: Vertical phase mixing: one-armed phase-spiral corresponding to
$n=1$ bending mode excited by the encounter with Sgr for MW disk+halo and MW
disk models (rows) at $R_{\rm c}=8\>{\rm kpc}$ and $12\>{\rm kpc}$ (columns).
The presence of DM halo slows down the rate of phase mixing, leading to more
loosely wrapped phase-spirals. Phase mixing occurs more rapidly in the inner
disk than in the outer disk.

The left panel of Fig. 9 plots the vertical frequency, $\Omega_{z}$, as a
function of the logarithm of the action, $I_{z}$, for the MW potential with
and without the halo and at guiding radii, $R_{\rm c}=8$ (red) and $12\>{\rm
kpc}$ (blue). The middle panel shows the behaviour of the corresponding phase
mixing timescale, $\tau_{\phi}$, as a function of $I_{z}$. Fig. 10 shows the
$(n,l,m)=(1,0,0)$ phase-spirals $400\>{\rm Myr}$ after the penultimate disk
crossing of Sagittarius, color coded by the MW disk response, $f_{1,100}$,
with blue (red) indicating higher (lower) phase-space density. Results for the
cases with and without the halo are shown in rows and for $R_{c}=8$ and
$12\>{\rm kpc}$ are shown in columns. Finally, the right panel of Fig. 9 shows
the $w_{z}=0$ cuts of the normalized response, $f_{1,100}/f_{0}$, as a
function of $I_{z}$, for the four different phase-spirals shown in Fig. 10.
The vertical frequency, $\Omega_{z}$, is a decreasing function of $\ln{I_{z}}$
in all cases, indicating that stars with larger actions (i.e., larger vertical
excursion amplitudes) oscillate slower. Note that $\left|{\rm
d}\Omega_{z}/{\rm d}\ln{I_{z}}\right|$ is an increasing (decreasing) function
of $I_{z}$ at small (large) $I_{z}$, reaching a maximum at intermediate
$I_{z}$. Consequently, the phase mixing timescale, $\tau_{\phi}$, which is the
inverse of $\left|{\rm d}\Omega_{z}/{\rm d}\ln{I_{z}}\right|$, attains its
minimum at $I_{z}/(h_{z}\sigma_{z})\sim 1$. Thus phase mixing occurs the
fastest at intermediate actions and slows down at larger actions, causing the
spiral to become more loosely wound (larger pitch angle) farther away from its
origin.

The rate of phase mixing is different in the four different cases. Closer to
the galactic center where the potential is deeper and steeper, stars have a
larger range of $\Omega_{z}$, or in other words $\Omega_{z}$ falls off more
steeply with $\ln{I_{z}}$ in the inner disk than in the outskirts. This leads
to faster phase mixing and therefore a much more tightly wound phase-spiral in
the inner disk (left panels of Fig. 10) as opposed to the outer disk (right
panels). The difference in the phase mixing rates is also manifest in the
$w_{z}=0$ response shown in the right panel of Fig. 9; note the longer
oscillation wavelengths of the blue lines (outer disk) as opposed to the red
lines (inner disk). Hence, in agreement with expectations, the inner part of
the disk equilibrates much faster than the outer part.

The presence of a DM halo deepens the potential well and thus boosts the
oscillation frequencies. But the halo also steepens the potential such that
the range of frequencies is reduced, i.e., $\Omega_{z}$ falls off more mildly
with $\ln{I_{z}}$ than in the disk only case. This leads to slower phase
mixing and therefore more loosely wound phase-spirals in the presence of the
halo (upper panels of Fig. 10) than in its absence (lower panels), the effect
being more pronounced in the outer (right panels) than in the inner (left
panels) disk. Equivalently, the $w_{z}=0$ response in the right panel of Fig.
9 shows longer wavelength wiggles in presence of the halo.

The above sensitivity of the phase mixing timescale to the detailed galaxy
potential implies that one can use phase-spirals to constrain it. One can
unwind the observed phase-spiral by adopting a form for the galactic
potential. Only for the correct potential will the spiral be properly unwound,
i.e., the pitch-angle, $\phi_{\rm I}$, go to zero for all $I_{z}$ (modulo
measurement errors) at the same time, $t_{0}$, in the past. This $t_{0}$ then
corresponds to the time elapsed since the maximum strength of the perturbation
that triggered the phase-spiral. However, this method to constrain the total
potential (disk plus dark matter) of the MW relies on the assumption of a
single, impulsive perturbation as the trigger. In reality, the phase-spirals
may have been impacted by multiple, overlapping perturbations and/or by large-
scale temporal fluctuations in the overall potential, which would severely
hamper this technique (Tremaine et al., 2022). We intend to investigate the
promise of phase-spirals as probes of the galactic potential for different
kinds of perturbation in future work.

##  6 Conclusion

In this paper, we have developed a linear perturbative formalism to analyze
the response of a realistic disk galaxy (characterized by a pseudo-isothermal
DF) embedded in an ambient spherical DM halo (modelled by an NFW profile) to
perturbations of diverse spatiotemporal nature: bars, spiral arms, and
encounters with satellite galaxies. Adopting the radial epicyclic
approximation, we perturb the FPE up to linear order (in action-angle space)
in presence of a perturbing potential, $\Phi_{\rm P}$, to compute the post-
perturbation linear response in the DF, $f_{1}$. Without self-gravity to
reinforce the response, the oscillations in the response phase mix away due to
an intrinsic spread in the frequencies of stars, giving rise to spiral
features in the phase-space distribution known as phase-spirals. Depending on
the timescale of $\Phi_{\rm P}$, different modes of disk oscillation,
corresponding to different phase-spiral structures, are excited. We summarize
our conclusions as follows:

  * •

Following an impulsive perturbation, the $(n,l,m)$ mode of the disk response
consists of stars oscillating with frequencies, $n\Omega_{z}$,
$l\Omega_{r}\approx l\kappa$ and $m\Omega_{\phi}$, along vertical, radial and
azimuthal directions respectively. Since the frequencies depend on the
actions, primarily on the vertical action $I_{z}$ and the angular momentum
$L_{z}$, the response phase mixes away, spawning phase-spirals. The dominant
modes of vertical oscillation are the anti-symmetric bending ($n=1$) and
symmetric breathing ($n=2$) modes, which induce initial dipolar and
quadrupolar perturbations in the $z-v_{z}$ or
$I_{z}\cos{w_{z}}-I_{z}\sin{w_{z}}$ phase-space. Over time these features are
phase-wrapped into one- and two-armed phase-spirals, respectively, due to the
variation of $\Omega_{z}$ with $I_{z}$.

  * •

Since $\Omega_{z}$ and $\Omega_{\phi}$ both depend on $L_{z}$, the amplitude
of the $I_{z}\cos{w_{z}}-I_{z}\sin{w_{z}}$ phase-spiral damps away over time,
typically as $\sim 1/t$ (equation 28), at a coarse-grained level, i.e., upon
marginalization over $L_{z}$. Therefore, in a realistic disk with ordered
motion, lateral mixing causes phase-spirals to damp out much slower than in
the isothermal slab with unconstrained lateral velocities discussed in Paper
I, where it occurs like a Gaussian in time.

  * •

Collisional diffusion due to scatterings of stars by GMCs, DM substructure,
etc. damps away the disk response to a perturbation, and therefore the phase-
spiral amplitude, at a fine-grained level. Typically, the diffusion in actions
is much more efficient than that in angles. The action gradients of the
response, which predominantly arise from the action dependence of the
oscillation frequencies, are erased by collisional diffusion, causing a super-
exponential damping of the response over a timescale, $\tau_{\rm D}^{({\rm
I})}$, which is $\sim 0.6-0.7\>{\rm Gyr}$ in the Solar neighborhood. The
diffusion timescale is shorter in the inner disk, for stars with smaller
$I_{z}$, and for higher-$n$ modes.

  * •

The response to a bar or spiral arm with a fixed pattern speed, $\Omega_{\rm
P}$, is dominated by the near-resonant stars ($\Omega_{\rm
res}=n\Omega_{z}+l\kappa+m(\Omega_{\phi}-\Omega_{\rm P})\approx 0$),
especially in the adiabatic regime (slowly evolving perturber amplitude).
Moreover, phase-mixing occurs gradually in the near-resonant parts of phase-
space. Most of the strong resonances are confined to the disk-plane, such as
the co-rotation ($n=l=0$) and Lindblad ($n=0,l=\pm 1,m=\pm 2$) resonances. For
a transient bar or spiral arm whose amplitude varies over time as
$\sim\exp{\left[-\omega^{2}_{0}t^{2}\right]}$, the response is maximal when
$\omega_{0}\sim\Omega_{\rm res}$. In the impulsive limit
($\omega_{0}\gg\Omega_{\rm res}$), the response is power-law suppressed, while
in the adiabatic limit ($\omega_{0}\ll\Omega_{\rm res}$) it is suppressed
(super)-exponentially.

  * •

For a thin disk, since $\Omega_{z}$ is very different from $\Omega_{\phi}$ and
$\kappa$, the vertical modes ($n\neq 0$) are generally not resonant with the
radial and azimuthal ones and thus undergo phase mixing. The strength of a
vertical mode primarily depends on the nature of the perturbing potential,
most importantly its timescale. Slower pulses trigger mainly bending $(n=1)$
modes, while faster pulses excite more pronounced breathing $(n=2)$ modes.
Therefore, a transient bar or spiral arm with amplitude
$\sim\exp{\left[-\omega^{2}_{0}t^{2}\right]}$ triggers a bending (breathing)
mode when the pulse-frequency, $\omega_{0}$, is smaller (larger) than
$\Omega_{z}$. The response to very slow perturbations
($\omega_{0}\ll\Omega_{z}$) is however heavily suppressed (adiabatic
shielding).

  * •

For a persistent bar or spiral arm with a fixed pattern speed, $\Omega_{\rm
P}$, that grows and saturates over time, the response initially develops a
phase-spiral. However, this transient response is quickly taken over by
coherent oscillations at the driving frequency, $\Omega_{\rm P}$, which
manifest in phase-space as a steadily rotating dipole (quadrupole) for the
bending (breathing) mode. Therefore, a transient (pulse-like) perturbation,
such as a bar or spiral arm whose amplitude varies over a timescale comparable
to the vertical oscillation period, $T_{z}\sim h_{z}/\sigma_{z}$, is essential
for the formation of a phase-spiral in $z$-$v_{z}$ space.

  * •

The above analysis suggests that if the recently discovered two-armed Gaia
phase-spiral (breathing mode) in the inner disk of the MW was indeed induced
by a spiral arm/bar as suggested by Hunt et al. (2022) using N-body
simulations, the spiral arm/bar was probably a transient one with a
predominantly symmetric vertical profile whose amplitude varied over a
timescale comparable to the vertical oscillation period. However, it remains
to be seen whether such a rapid excitation and decay of a spiral arm/bar
perturbation is realistic.

  * •

We have computed the response of the MW disk, embedded in an extended DM halo,
to disk-crossing perturbations by several of its satellite galaxies. We find
that the response in the Solar neighborhood is dominated by the perturbations
due to Sgr, followed by those due to the LMC, Hercules and Leo II. This
implies that, if the Gaia snail near the Solar radius was indeed triggered by
a MW satellite (which is still subject to debate), Sgr is the leading
contender (see also Banik et al., 2022). However, if that is the case, then
the impact (disk crossing) must have happened within the last $\sim
0.6-0.7\>{\rm Gyr}$ in order for the response to have survived damping due to
collisional diffusion.

  * •

The amplitude of the response (at a fixed guiding radius $R_{\rm c}$) to
satellite encounters scales as $v_{\rm P}^{-1}$ in the impulsive (large
$v_{\rm P}$) limit, but is exponentially suppressed in the adiabatic (small
$v_{\rm P}$) limit, a phenomenon known as adiabatic shielding. The resonant
modes with $n\Omega_{z}+l\kappa+m\Omega_{\phi}=0$ are not suppressed but
rather become non-linear in the adiabatic regime. The peak response of a mode
(with $n\Omega_{z}+l\kappa+m\Omega_{\phi}\neq 0$) is achieved at intermediate
velocities for which the encounter frequency is commensurate with the
oscillation frequencies of the stars, i.e., the near-resonance condition given
by equation (42) is satisfied.

  * •

The response of a disk to an encounter with a satellite galaxy depends
primarily on three parameters: (i) impact velocity $v_{\rm P}$, (ii) polar
angle of impact $\theta_{\rm P}$, and (iii) position on the disk relative to
the point of impact where the satellite crossed the disk. Slower (faster)
encounters excite predominantly $n=1$ bending ($n=2$ breathing) modes. More
planar encounters (those with larger $\theta_{\rm P}$) typically result in
larger breathing-to-bending ratios farther away from the impact radius while
this trend gets reversed closer to it. In general, breathing modes dominate
over bending modes closer to the point of impact, in agreement with $N$-body
simulations of the MW-Sgr encounter (Hunt et al., 2021). Since the impact
velocities of the MW satellites are all fairly similar to the local circular
velocity, the decisive factor for breathing vs. bending modes is not so much
the impact velocity, but rather the distance from the point of impact.

  * •

The $(n,m)=(1,-2)$ and $(-1,2)$ modes generally dominate the response for
slower satellite encounters, e.g., that of Sgr with respect to the Solar
neighborhood, due to the tidal distortion of the disk by the satellite. The
in-plane spatial structure of the disk response therefore generally resembles
a two-armed warp ($n=1$) or spiral ($n=2$).

  * •

The presence of an extended DM halo causes phase mixing to occur slower, and
modifies the structural appearance of the phase-spirals (i.e., the pitch angle
as function of vertical action). Hence, provided that the phase-spiral was
triggered by a single, impulsive perturbation, the detailed shape of the
phase-spiral can in principle be used to constrain not only the time elapsed
since the perturbation (Darragh-Ford et al., 2023) but also the total
(disk$+$halo) potential. If the phase-spiral has been triggered and/or
impacted by multiple, overlapping perturbations the situation is less clear.
In future work we intend to investigate the constraining power of phase-
spirals for different kinds of perturbation.

This paper focused on the analysis of the phase mixing component (phase-
spirals) of the ‘direct’ disk response to various perturbations such as bars,
spiral arms and satellite galaxies. However this leaves out some other
potentially important features of the disk response. First of all, we
considered the ambient DM halo to be non-responsive. In reality, the DM halo
will also be perturbed, for example by an impacting satellite, and this halo
response, which can be enhanced by self-gravity, can indirectly perturb the
disk. A preliminary, perturbative analysis based on the $N$-body simulation of
the MW-Sgr encounter by Hunt et al. (2021) suggests that the indirect disk
response to halo perturbations (triggered by Sgr) is comparable, but sub-
dominant, to the direct response to Sgr. However, a more detailed analysis is
warranted, which we leave for future work. Secondly, we have neglected the
self-gravity of the disk response. As discussed in Paper I, the dominant
effect of self-gravity is to cause coherent point mode oscillations (Mathur,
1990; Weinberg, 1991; Darling & Widrow, 2019b) of the disk, which in linear
theory are decoupled from the phase-spirals. However, self-gravity can enhance
the amplitude of the phase-spiral. Although recent developments (Dootson &
Magorrian, 2022) have shed some light on the self-gravitating response of
razor-thin disks to bar perturbations, a more generic theoretical description
of the self-gravitating response of inhomogeneous, thick disks to general
perturbations (bars, spiral arms, satellite galaxies, etc) is still lacking.
We hope to include the effects of self-gravity on disk perturbations in future
work.

## Acknowledgments

The authors are grateful to Kathryn Johnston, Jason Hunt, Adrian Price-Whelan,
Chris Hamilton, Elise Darragh-Ford, James Binney, John Magorrian and Elena
D’Onghia for insightful discussions and valuable suggestions. FCvdB is
supported by the National Aeronautics and Space Administration through Grant
No. 19-ATP19-0059 issued as part of the Astrophysics Theory Program. MDW is
supported by the National Science Foundation through Grant No. AST-1812689.

## References

  * Antoja et al. (2018) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360, doi: 10.1038/s41586-018-0510-7
  * Banik & van den Bosch (2021a) Banik, U., & van den Bosch, F. C. 2021a, ApJ, 912, 43, doi: 10.3847/1538-4357/abeb6d
  * Banik & van den Bosch (2021b) —. 2021b, MNRAS, 502, 1441, doi: 10.1093/mnras/stab092
  * Banik & van den Bosch (2022) —. 2022, ApJ, 926, 215, doi: 10.3847/1538-4357/ac4242
  * Banik et al. (2022) Banik, U., Weinberg, M. D., & van den Bosch, F. C. 2022, ApJ, 935, 135, doi: 10.3847/1538-4357/ac7ff9
  * Bekki & Stanimirović (2009) Bekki, K., & Stanimirović, S. 2009, MNRAS, 395, 342, doi: 10.1111/j.1365-2966.2009.14514.x
  * Bennett et al. (2022) Bennett, M., Bovy, J., & Hunt, J. A. S. 2022, ApJ, 927, 131, doi: 10.3847/1538-4357/ac5021
  * Binney (2010) Binney, J. 2010, MNRAS, 401, 2318, doi: 10.1111/j.1365-2966.2009.15845.x
  * Binney & Lacey (1988) Binney, J., & Lacey, C. 1988, MNRAS, 230, 597, doi: 10.1093/mnras/230.4.597
  * Binney & Schönrich (2018) Binney, J., & Schönrich, R. 2018, MNRAS, 481, 1501, doi: 10.1093/mnras/sty2378
  * Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton University Press)
  * Binney & Tremaine (2008) —. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bland-Hawthorn & Tepper-García (2021) Bland-Hawthorn, J., & Tepper-García, T. 2021, MNRAS, 504, 3168, doi: 10.1093/mnras/stab704
  * Bland-Hawthorn et al. (2019) Bland-Hawthorn, J., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 486, 1167, doi: 10.1093/mnras/stz217
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Bovy & Rix (2013) Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115, doi: 10.1088/0004-637X/779/2/115
  * Carlberg & Sellwood (1985) Carlberg, R. G., & Sellwood, J. A. 1985, ApJ, 292, 79, doi: 10.1086/163134
  * Chiba & Schönrich (2022) Chiba, R., & Schönrich, R. 2022, MNRAS, 513, 768, doi: 10.1093/mnras/stac697
  * Darling & Widrow (2019a) Darling, K., & Widrow, L. M. 2019a, MNRAS, 490, 114, doi: 10.1093/mnras/stz2539
  * Darling & Widrow (2019b) —. 2019b, MNRAS, 484, 1050, doi: 10.1093/mnras/sty3508
  * Darragh-Ford et al. (2023) Darragh-Ford, E., Hunt, J. A. S., Price-Whelan, A. M., & Johnston, K. V. 2023, arXiv e-prints, arXiv:2302.09086, doi: 10.48550/arXiv.2302.09086
  * Dootson & Magorrian (2022) Dootson, D., & Magorrian, J. 2022, arXiv e-prints, arXiv:2205.15725.  https://arxiv.org/abs/2205.15725
  * Erkal et al. (2019) Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685, doi: 10.1093/mnras/stz1371
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gandhi et al. (2022) Gandhi, S. S., Johnston, K. V., Hunt, J. A. S., et al. 2022, ApJ, 928, 80, doi: 10.3847/1538-4357/ac47f7
  * Gnedin & Ostriker (1999) Gnedin, O. Y., & Ostriker, J. P. 1999, ApJ, 513, 626, doi: 10.1086/306864
  * Gradshteyn & Ryzhik (1965) Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of integrals, series and products
  * Hamilton et al. (2022) Hamilton, C., Tolman, E. A., Arzamasskiy, L., & Duarte, V. N. 2022, arXiv e-prints, arXiv:2208.03855.  https://arxiv.org/abs/2208.03855
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hunt et al. (2022) Hunt, J. A. S., Price-Whelan, A. M., Johnston, K. V., & Darragh-Ford, E. 2022, MNRAS, doi: 10.1093/mnrasl/slac082
  * Hunt et al. (2021) Hunt, J. A. S., Stelea, I. A., Johnston, K. V., et al. 2021, MNRAS, 508, 1459, doi: 10.1093/mnras/stab2580
  * Johnston et al. (2017) Johnston, K. V., Price-Whelan, A. M., Bergemann, M., et al. 2017, Galaxies, 5, 44, doi: 10.3390/galaxies5030044
  * Kaur & Sridhar (2018) Kaur, K., & Sridhar, S. 2018, ApJ, 868, 134, doi: 10.3847/1538-4357/aaeacf
  * Kaur & Stone (2022) Kaur, K., & Stone, N. C. 2022, MNRAS, 515, 407, doi: 10.1093/mnras/stac1729
  * Laporte et al. (2018) Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286, doi: 10.1093/mnras/sty1574
  * Laporte et al. (2019) Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134, doi: 10.1093/mnras/stz583
  * Li & Widrow (2021) Li, H., & Widrow, L. M. 2021, MNRAS, 503, 1586, doi: 10.1093/mnras/stab574
  * Li (2021) Li, Z.-Y. 2021, ApJ, 911, 107, doi: 10.3847/1538-4357/abea17
  * Li et al. (2020) Li, Z.-Z., Qian, Y.-Z., Han, J., et al. 2020, ApJ, 894, 10, doi: 10.3847/1538-4357/ab84f0
  * Lichtenberg & Lieberman (1992) Lichtenberg, A., & Lieberman, M. 1992, Regular and Chaotic Dynamics
  * Łokas (2009) Łokas, E. L. 2009, MNRAS, 394, L102, doi: 10.1111/j.1745-3933.2009.00620.x
  * Lynden-Bell (1962) Lynden-Bell, D. 1962, MNRAS, 124, 279, doi: 10.1093/mnras/124.4.279
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1, doi: 10.1093/mnras/157.1.1
  * Maoz (1991) Maoz, E. 1991, ApJ, 375, 687, doi: 10.1086/170232
  * Mathur (1990) Mathur, S. D. 1990, MNRAS, 243, 529
  * McMillan (2011) McMillan, P. J. 2011, MNRAS, 414, 2446, doi: 10.1111/j.1365-2966.2011.18564.x
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Price-Whelan et al. (2020) Price-Whelan, A., Sipőcz, B., Lenz, D., et al. 2020, adrn/gala: v1.3, v1.3, Zenodo, doi: 10.5281/zenodo.4159870
  * Price-Whelan (2017) Price-Whelan, A. M. 2017, The Journal of Open Source Software, 2, doi: 10.21105/joss.00388
  * Riley et al. (2019) Riley, A. H., Fattahi, A., Pace, A. B., et al. 2019, MNRAS, 486, 2679, doi: 10.1093/mnras/stz973
  * Simon & Geha (2007) Simon, J. D., & Geha, M. 2007, ApJ, 670, 313, doi: 10.1086/521816
  * Smith et al. (2015) Smith, R., Flynn, C., Candlish, G. N., Fellhauer, M., & Gibson, B. K. 2015, MNRAS, 448, 2934, doi: 10.1093/mnras/stv228
  * Sridhar (1989) Sridhar, S. 1989, MNRAS, 238, 1159, doi: 10.1093/mnras/238.4.1159
  * Tremaine et al. (2022) Tremaine, S., Frankel, N., & Bovy, J. 2022, arXiv e-prints, arXiv:2212.11990.  https://arxiv.org/abs/2212.11990
  * Tremaine & Weinberg (1984) Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729, doi: 10.1093/mnras/209.4.729
  * van den Bosch et al. (2018) van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043, doi: 10.1093/mnras/stx2956
  * Vasiliev & Belokurov (2020) Vasiliev, E., & Belokurov, V. 2020, MNRAS, 497, 4162, doi: 10.1093/mnras/staa2114
  * Weinberg (1989) Weinberg, M. D. 1989, MNRAS, 239, 549, doi: 10.1093/mnras/239.2.549
  * Weinberg (1991) —. 1991, ApJ, 373, 391, doi: 10.1086/170059
  * Weinberg (1994a) —. 1994a, AJ, 108, 1398, doi: 10.1086/117161
  * Weinberg (1994b) —. 1994b, AJ, 108, 1403, doi: 10.1086/117162
  * Weinberg (2004) —. 2004, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0404169
  * Widmark et al. (2022a) Widmark, A., Hunt, J. A. S., Laporte, C. F. P., & Monari, G. 2022a, A&A, 663, A16, doi: 10.1051/0004-6361/202243173
  * Widmark et al. (2022b) Widmark, A., Laporte, C. F. P., & Monari, G. 2022b, A&A, 663, A15, doi: 10.1051/0004-6361/202142819
  * Widrow et al. (2014) Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971, doi: 10.1093/mnras/stu396

##  Appendix A The unperturbed galaxy

Under the radial epicyclic approximation (small $I_{R}$), the unperturbed DF,
$f_{0}$, for a rotating MW-like disk galaxy can be well approximated as a
pseudo-isothermal DF, i.e., written as a nearly isothermal separable function
of the azimuthal, radial and vertical actions. Following Binney (2010), we
write

$\displaystyle
f_{0}=\frac{1}{\pi}{\left(\frac{\Omega_{\phi}\Sigma}{\kappa\,\sigma^{2}_{R}}\right)}_{R_{\rm
c}}\left(1+\tanh{\frac{L_{z}}{L_{0}}}\right)\times\exp{\left[-\frac{\kappa
I_{R}}{\sigma^{2}_{R}}\right]}\times\frac{1}{\sqrt{2\pi}h_{z}\sigma_{z}}\exp{\left[-\frac{E_{z}(I_{z})}{\sigma^{2}_{z}}\right]}\,.$
(A1)

The vertical structure of this disk is isothermal, while the radial profile is
pseudo-isothermal. Here $\Sigma=\Sigma(R)$ is the surface density of the disk,
$L_{z}$ is the $z$-component of the angular momentum, which is equal to
$I_{\phi}$, $R_{\rm c}=R_{\rm c}(L_{z})$ is the guiding radius,
$\Omega_{\phi}$ is the circular frequency, and $\kappa=\kappa(R_{\rm
c})=\lim_{I_{R}\to 0}{\Omega_{R}}$ is the radial epicyclic frequency (Binney &
Tremaine, 1987). If $L_{0}$ is sufficiently small, then we can further
approximate the above form for $f_{0}$ as

$\displaystyle
f_{0}\approx\frac{\sqrt{2}}{\pi^{3/2}\,\sigma_{z}h_{z}}{\left(\frac{\Omega_{\phi}\Sigma}{\kappa\,\sigma^{2}_{R}}\right)}_{R_{\rm
c}}\,\exp{\left[-\frac{\kappa
I_{R}}{\sigma^{2}_{R}}\right]}\,\exp{\left[-\frac{E_{z}(I_{z})}{\sigma^{2}_{z}}\right]}\,\Theta(L_{z})\,,$
(A2)

where $\Theta(x)$ is the Heaviside step function. Thus we assume that the
entire galaxy is composed of prograde stars with $L_{z}>0$.

The corresponding density profile can be written as a product of an
exponential radial profile and an isothermal ($\operatorname{sech}^{2}$)
vertical profile, i.e.,

$\displaystyle\rho(R,z)=\rho_{\rm
c}\,\exp{\left[-\frac{R}{h_{R}}\right]}\,\operatorname{sech}^{2}{\left(\frac{z}{h_{z}}\right)},$
(A3)

where $h_{R}$ and $h_{z}$ are the radial and vertical scale heights,
respectively. Throughout we adopt the thin disk limit, i.e., $h_{z}\ll h_{R}$.
The surface density profile is given by

$\displaystyle\Sigma(R)=\int_{-\infty}^{\infty}{\rm d}z\,\rho(R,z)=\Sigma_{\rm
c}\exp{\left[-\frac{R}{h_{R}}\right]},$ (A4)

where $\Sigma_{\rm c}=\rho_{\rm c}h_{z}$ is the central surface density of the
disk. We assume a radially varying vertical velocity dispersion, $\sigma_{z}$,
satisfying $\sigma^{2}_{z}(R)=2\pi Gh_{z}\Sigma(R)$ (Binney & Tremaine, 2008).
We assume a similar profile for $\sigma^{2}_{R}$ such that the ratio,
$\sigma_{R}/\sigma_{z}$ is constant throughout the disk (Binney, 2010) and
equal to the value at the Solar vicinity.

Throughout this paper, for the ease of computation of the frequencies (because
of a simple analytic form of the potential), we approximate the above density
profile by a combination of three Miyamoto & Nagai (1975) disk profiles (Smith
et al., 2015), i.e., the 3MN profile as implemented in the Gala Python package
(Price-Whelan, 2017; Price-Whelan et al., 2020). The corresponding disk
potential is given by

$\displaystyle\Phi_{\rm
d}(R,z)=-\sum_{i=1}^{3}\frac{GM_{i}}{\sqrt{R^{2}+{\left(a_{i}+\sqrt{z^{2}+b^{2}_{i}}\right)}^{2}}},$
(A5)

where $M_{i}$, $a_{i}$ and $b_{i}$, with $i=1,2,3$, are the mass, scale radius
and scale height corresponding to each of the MN profiles.

The MW disk is believed to be embedded in a much more extended DM halo, which
we model using a spherical NFW (Navarro et al., 1997) profile with potential

$\displaystyle\Phi_{\rm h}(R,z)=-\frac{GM_{\rm vir}}{R_{\rm
vir}}\,\frac{c}{f(c)}\,\frac{\ln(1+r/r_{\rm s})}{r/r_{\rm s}}.$ (A6)

Here $M_{\rm vir}$ is the virial mass of the halo, $r_{s}$ is the scale
radius, $c=R_{\rm vir}/r_{s}$ is the concentration ($R_{\rm vir}$ is the
virial radius), and $f(c)=\ln{\left(1+c\right)}-c/(1+c)$. The combined
potential experienced by the disk stars is thus given by

$\displaystyle\Phi_{0}(R,z)=\Phi_{\rm d}(R,z)+\Phi_{\rm h}(R,z).$ (A7)

##  Appendix B Fourier coefficients of spiral arm or bar perturbing potential

An essential ingredient of the disk response to spiral arm or bar
perturbations is the Fourier component of the perturber potential,
$\Phi_{nlm}$. This can be computed as follows:

$\displaystyle\Phi_{nlm}({\bf I},t)$
$\displaystyle=\frac{1}{{\left(2\pi\right)}^{3}}\int_{0}^{2\pi}{\rm
d}w_{z}\int_{0}^{2\pi}{\rm d}w_{R}\int_{0}^{2\pi}{\rm
d}w_{\phi}\,\exp{\left[-i(nw_{z}+lw_{R}+mw_{\phi})\right]}\,\Phi_{\rm
P}\left({\bf r},t\right).$ (B1)

To evaluate this first we need to calculate ${\bf r}=(z,R,\phi)$ as a function
of $\left({\bf w},{\bf I}\right)=(w_{z},w_{\phi},w_{R},{\bf I})$ where
$I_{\phi}=L_{z}$, the angular momentum. Under the epicyclic approximation, $R$
can be expressed as a sum of the guiding radius and an oscillating epicyclic
term, i.e.,

$\displaystyle R\approx R_{\rm
c}(L_{z})+\sqrt{\frac{2I_{R}}{\kappa}}\sin{w_{R}},$ (B2)

and the azimuthal angle, $w_{\phi}$, is given by

$\displaystyle w_{\phi}\approx\phi-\frac{2\,\Omega_{\phi}}{R_{\rm
c}\kappa}\sqrt{\frac{2I_{R}}{\kappa}}\cos{\theta_{R}}.$ (B3)

The vertical distance $z$ from the mid-plane is related to $R_{\rm c}(L_{z})$
and $\left(w_{z},I_{z}\right)$, according to

$\displaystyle w_{z}=\Omega_{z}(R_{\rm c},I_{z})\int_{0}^{z}\frac{{\rm
d}z^{\prime}}{\sqrt{2\left[E_{z}(R_{\rm c},I_{z})-\Phi_{z}(R_{\rm
c},z^{\prime})\right]}},$ (B4)

where $\Omega_{z}(R_{\rm c},I_{z})=2\pi/T_{z}(R_{\rm c},I_{z})$, with
$T_{z}(R_{\rm c},I_{z})$ given by Equation (19). The above equation can be
numerically inverted to obtain $z(R_{\rm c},w_{z},I_{z})$.

Upon substituting the above expressions for $R$, $\phi$ and $z$ in terms of
$({\bf w},{\bf I})$ in the expression for $\Phi_{\rm P}$ given in equation
(21), we obtain

$\displaystyle\Phi_{nlm}\left({\bf I},t\right)$ $\displaystyle=-\frac{2\pi
G\Sigma_{\rm
P}}{k_{R}}\left(\sum_{m_{\phi}=0,2,-2}\delta_{m,m_{\phi}}\right)\frac{\operatorname{sgn}(m)\exp{\left[i\,{\rm
sgn}(m)k_{R}R_{\rm c}(I_{\phi})\right]}}{2i}$
$\displaystyle\times\exp{\left[i\,l\tan^{-1}{\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\sqrt{\frac{2I_{R}}{\kappa}}}\right]}\,J_{l}\left(\sqrt{k_{R}^{2}+{\left(\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\right)}^{2}}\sqrt{\frac{2I_{R}}{\kappa}}\right)$
$\displaystyle\times\left[\alpha\,{\cal M}_{\rm o}(t)\Phi_{n}^{({\rm
o})}(I_{z})+{\cal M}_{\rm e}(t)\Phi_{n}^{({\rm
e})}(I_{z})\right]\exp{\left[-im\Omega_{\rm P}t\right]},$ (B5)

where $J_{l}$ is the $l^{\rm th}$ order Bessel function of the first kind,

$\displaystyle\operatorname{sgn}(m)$ $\displaystyle=\begin{cases}1,&m\geq
0,\\\ -1,&m<0,\end{cases}$ (B6)

and $\Phi_{n}^{({\rm o})}(I_{z})$ and $\Phi_{n}^{({\rm e})}(I_{z})$ are given
by

$\displaystyle\Phi_{n}^{({\rm o})}(I_{z})$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm d}w_{z}\sin{nw_{z}}\,{\cal
F}_{\rm o}\left(z,k_{z}^{({\rm o})}\right),$ $\displaystyle\Phi_{n}^{({\rm
e})}(I_{z})$ $\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}{\rm
d}w_{z}\cos{nw_{z}}\,{\cal F}_{\rm e}\left(z,k_{z}^{({\rm e})}\right).$ (B7)

In deriving equation (B5) we have used the Hansen-Bessel formula which
provides the following integral representation for Bessel functions of the
first kind,

$\displaystyle\int_{0}^{2\pi}{\rm
d}x\,\exp{\left[-ilx\right]}\,\exp{\left[i\alpha\sin{x}\right]}=2\pi
J_{l}\left(\alpha\right),$ (B8)

and the identity for expansion in products of Bessel functions given in
equation (8.530.2) of Gradshteyn & Ryzhik (1965, see also section 6.1 of ()).
We have also used the identity,

$\displaystyle\int_{0}^{2\pi}{\rm
d}\phi\,\exp{\left[-im\phi\right]}=2\pi\,\delta_{m,0}.$ (B9)

##  Appendix C Perturbation by encounter with satellite galaxy

###  C.1 Computation of the disk response

To evaluate the disk response to satellite encounters using equation (20) we
first evaluate the $\tau$ integral (with $t_{\rm i}\to-\infty$) of the
satellite potential given in equation (37) and then compute the Fourier
transform of the result. This yields the expression for the response in
equation (20) with

$\displaystyle{\cal I}_{nlm}({\bf I},t)=\exp{\left[-i\Omega
t\right]}\int_{-\infty}^{t}{\rm
d}\tau\,\exp{\left[i\Omega\tau\right]}\,\Phi_{nlm}\left({\bf I},\tau\right)$
$\displaystyle=\frac{\exp{\left[-i\Omega
t\right]}}{{\left(2\pi\right)}^{3}}\int_{0}^{2\pi}{\rm
d}w_{z}\exp{\left[-inw_{z}\right]}\int_{0}^{2\pi}{\rm
d}w_{R}\exp{\left[-ilw_{R}\right]}\int_{0}^{2\pi}{\rm
d}w_{\phi}\exp{\left[-imw_{\phi}\right]}\int_{-\infty}^{t}{\rm
d}\tau\,\exp{\left[i\Omega\tau\right]}\,\Phi_{\rm P}(z,R,\phi,\tau),$ (C1)

where

$\displaystyle\Omega=n\Omega_{z}+l\Omega_{R}+m\Omega_{\phi}.$ (C2)

We perform the inner $\tau$ integral of $\Phi_{\rm P}$ to obtain

$\displaystyle\int_{-\infty}^{t}{\rm
d}\tau\,\exp{\left[i\Omega\tau\right]}\,\Phi_{\rm P}(z,R,\phi,\tau)$
$\displaystyle=-\frac{GM_{\rm P}}{v_{\rm P}}\exp{\left[i\frac{\Omega{\cal
S}}{v_{\rm P}}\right]}\int_{-\infty}^{t-{\cal S}/v_{\rm P}}{\rm
d}\tau\frac{\exp{\left[i\Omega\tau\right]}}{\sqrt{\tau^{2}+{\left({\cal
R}^{2}+\varepsilon^{2}\right)}/{v^{2}_{\rm P}}}}$
$\displaystyle=-\frac{GM_{\rm P}}{v_{\rm P}}\exp{\left[i\frac{\Omega{\cal
S}}{v_{\rm P}}\right]}\int_{-\infty}^{\left(v_{\rm P}t-{\cal
S}\right)/\sqrt{{\cal R}^{2}+\varepsilon^{2}}}{\rm
d}x\frac{\exp{\left[i\left(\Omega\sqrt{{\cal R}^{2}+\varepsilon^{2}}/v_{\rm
P}\right)x\right]}}{\sqrt{x^{2}+1}}$ $\displaystyle=-\frac{2\,GM_{\rm
P}}{v_{\rm P}}\exp{\left[i\frac{\Omega{\cal S}}{v_{\rm
P}}\right]}K_{0i}\left(\frac{\Omega\sqrt{{\cal R}^{2}+\varepsilon^{2}}}{v_{\rm
P}},\frac{v_{\rm P}t-{\cal S}}{\sqrt{{\cal R}^{2}+\varepsilon^{2}}}\right).$
(C3)

Here $K_{0i}$ is defined as

$\displaystyle K_{0i}(\alpha,\beta)=\frac{1}{2}\int_{-\infty}^{\beta}{\rm
d}x\,\frac{\exp{\left[i\alpha x\right]}}{\sqrt{x^{2}+1}},$ (C4)

which asymptotes to the zero-th order modified Bessel function of the second
kind, $K_{0}\left(\left|\alpha\right|\right)$, in the limit $\beta\to\infty$.
${\cal R}$ and ${\cal S}$ are respectively the projections perpendicular and
parallel to the direction of ${\bf v_{\rm P}}$ of the vector connecting the
point of observation, $(z,R,\phi)$, with the point of impact, and are given by

$\displaystyle{\cal R}^{2}$ $\displaystyle={\left[R\sin{(\phi-\phi_{\rm
P})}+R_{\rm d}\sin{\phi_{\rm P}}\right]}^{2}+{\left[(R\cos{(\phi-\phi_{\rm
P})}-R_{\rm d}\cos{\phi_{\rm P}})\cos{\theta_{\rm P}}-z\sin{\theta_{\rm
P}}\right]}^{2}$ $\displaystyle{\cal S}$ $\displaystyle=(R\cos{(\phi-\phi_{\rm
P})}-R_{\rm d}\cos{\phi_{\rm P}})\sin{\theta_{\rm P}}+z\cos{\theta_{\rm P}}.$
(C5)

In deriving equation (C3), we have only considered the direct term in the
expression for $\Phi_{\rm P}$ given in equation (37); the indirect term turns
out to be sub-dominant.

In the large time limit, i.e., $t\gg{\cal S}/v_{\rm P}$, $K_{0i}$ asymptotes
to $K_{0}\left(\left|\Omega\right|\sqrt{{\cal R}^{2}+\varepsilon^{2}}/v_{\rm
P}\right)$. We substitute the expressions for $R$ and $z$ in terms of $({\bf
w},{\bf I})$ given in equations (B2) and (B4) in the above expressions for
${\cal R}$ and ${\cal S}$. Further substituting the resultant $\tau$ integral
from equation (C3) in equation (C1), substituting $w_{\phi}$ in terms of
$\phi$ using equation (B3), adopting the small $I_{R}$ limit, and performing
the $w_{R}$ integral, we obtain

$\displaystyle{\cal I}_{nlm}({\bf I},t)\approx-\frac{2GM_{\rm P}}{v_{\rm
P}}\exp{\left[-i\Omega
t\right]}\times\exp{\left[-i\frac{\Omega\sin{\theta_{\rm P}}\cos{\phi_{\rm
P}}}{v_{\rm P}}R_{\rm
d}\right]}\times\exp{\left[i\,l\tan^{-1}{\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\sqrt{\frac{2I_{R}}{\kappa}}}\right]}$
$\displaystyle\times\frac{1}{{\left(2\pi\right)}^{2}}\int_{0}^{2\pi}dw_{z}\exp{\left[-inw_{z}\right]}\exp{\left[i\frac{\Omega\cos{\theta_{\rm
P}}}{v_{\rm
P}}z\right]}\int_{0}^{2\pi}d\phi\,\exp{\left[-im\phi\right]}\,\exp{\left[i\frac{\Omega\sin{\theta_{\rm
P}}\cos{\left(\phi-\phi_{\rm P}\right)}}{v_{\rm P}}R_{\rm c}\right]}$
$\displaystyle\times J_{l}\left(\sqrt{{\left(\frac{\Omega\sin{\theta_{\rm
P}}}{v_{\rm P}}\right)}^{2}\cos^{2}{\left(\phi-\phi_{\rm
P}\right)}+{\left(\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\right)}^{2}}\sqrt{\frac{2I_{R}}{\kappa}}\right)K_{0i}\left(\frac{\Omega\sqrt{{\cal
R}^{2}_{\rm c}+\varepsilon^{2}}}{v_{\rm P}},\frac{v_{\rm P}t-{\cal S}_{\rm
c}}{\sqrt{{\cal R}^{2}_{\rm c}+\varepsilon^{2}}}\right),$ (C6)

where ${\cal R}_{\rm c}={\cal R}(R=R_{\rm c})$ and ${\cal S}_{\rm c}={\cal
S}(R=R_{\rm c})$. Here we have used the integral representation of Bessel
functions of the first kind given in equation (B8) and the identity given in
equation (8.530.2) of Gradshteyn & Ryzhik (1965).

The expression for ${\cal I}_{nlm}$ given in equation (C6) consists of the
leading order expansion in $\sqrt{2I_{R}/\kappa}$. A more precise expression
that is accurate up to second order in $\sqrt{2I_{R}/\kappa}$ is given, in the
large time limit, as

$\displaystyle{\cal I}_{nlm}({\bf
I},t)\overset{\mathrm{t\to\infty}}{\approx}-\frac{2GM_{\rm P}}{v_{\rm
P}}\exp{\left[-i\Omega
t\right]}\times\exp{\left[-i\frac{\Omega\sin{\theta_{\rm P}}\cos{\phi_{\rm
P}}}{v_{\rm P}}R_{\rm d}\right]}$
$\displaystyle\times\frac{1}{{\left(2\pi\right)}^{2}}\int_{0}^{2\pi}dw_{z}\exp{\left[-inw_{z}\right]}\exp{\left[i\frac{\Omega\cos{\theta_{\rm
P}}}{v_{\rm
P}}z\right]}\int_{0}^{2\pi}d\phi\,\exp{\left[-im\phi\right]}\,\exp{\left[i\frac{\Omega\sin{\theta_{\rm
P}}\cos{\left(\phi-\phi_{\rm P}\right)}}{v_{\rm P}}R_{\rm c}\right]}$
$\displaystyle\times\exp{\left[i\,l\tan^{-1}{\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\sqrt{\frac{2I_{R}}{\kappa}}}\right]}\,\left[\zeta^{(0)}J_{l}\left(\chi\right)-i\zeta^{(1)}J^{\prime}_{l}\left(\chi\right)-\frac{1}{2}\zeta^{(2)}J^{\prime\prime}_{l}\left(\chi\right)\right],$
(C7)

where

$\displaystyle\chi=\sqrt{{\left(\frac{\Omega\sin{\theta_{\rm P}}}{v_{\rm
P}}\right)}^{2}\cos^{2}{\left(\phi-\phi_{\rm
P}\right)}+{\left(\frac{2m\Omega_{\phi}}{R_{\rm
c}\kappa}\right)}^{2}}\sqrt{\frac{2I_{R}}{\kappa}},$ (C8)

and

$\displaystyle\zeta^{(0)}$ $\displaystyle=K_{0}\left(\eta\right),$
$\displaystyle\zeta^{(1)}$
$\displaystyle=\sqrt{\frac{2I_{R}}{\kappa}}\,\frac{\partial{\cal R}_{\rm
c}}{\partial R_{\rm c}}\frac{{\cal R}_{\rm c}}{\sqrt{{\cal R}^{2}_{\rm
c}+\varepsilon^{2}}}\frac{\left|\Omega\right|}{v_{\rm
P}}K^{\prime}_{0}\left(\eta\right),$ $\displaystyle\zeta^{(2)}$
$\displaystyle=\frac{2I_{R}}{\kappa}\,\left[{\left(\frac{\partial{\cal R}_{\rm
c}}{\partial R_{\rm c}}\right)}^{2}\frac{{\cal R}^{2}_{\rm c}}{{\cal
R}^{2}_{\rm c}+\varepsilon^{2}}\frac{\Omega^{2}}{v^{2}_{\rm
P}}K^{\prime\prime}_{0}(\eta)+\left\\{\frac{\partial^{2}{\cal R}_{\rm
c}}{\partial R^{2}_{c}}\frac{{\cal R}_{\rm c}}{\sqrt{{\cal R}^{2}_{\rm
c}+\varepsilon^{2}}}+{\left(\frac{\partial{\cal R}_{\rm c}}{\partial R_{\rm
c}}\right)}^{2}\frac{\varepsilon^{2}}{{\left({\cal R}^{2}_{\rm
c}+\varepsilon^{2}\right)}^{3/2}}\right\\}\frac{\left|\Omega\right|}{v_{\rm
P}}K^{\prime}_{0}(\eta)\right],$ (C9)

with

$\displaystyle\eta$ $\displaystyle=\frac{\left|\Omega\right|\sqrt{{\cal
R}^{2}_{\rm c}+\varepsilon^{2}}}{v_{\rm P}}.$ (C10)

Here each prime denotes a single derivative of the function with respect to
its argument.

Substituting the expression for ${\cal I}_{nlm}$ given in equation (C7) in the
expression for the disk response given in equation (20), and adopting the
fiducial parameters for the MW galaxy and those corresponding to satellite
encounters (as detailed in section 4.2), we compute the response of the MW
disk to past and future encounters with its satellite galaxies. Results for
the steady state disk response (in the collisionless limit) of the
$(n,l,m)=(1,0,0)$ mode, corresponding to
$I_{z}=I_{z,\odot}=h_{z}\sigma_{z,\odot}$ and $R_{c}(L_{z})=8\>{\rm kpc}$ and
marginalized over $I_{R}$, are summarised in Table 1.

MW satellite | Mass | $f_{1,n=1}/f_{0}$ | $t_{\rm cross}$ | $f_{1,n=1}/f_{0}$ | $t_{\rm cross}$ | $f_{1,n=1}/f_{0}$ | $t_{\rm cross}$
---|---|---|---|---|---|---|---
name | $(\>{\rm M_{\odot}})$ |  | $(\>{\rm Gyr})$ |  | $(\>{\rm Gyr})$ |  | $(\>{\rm Gyr})$
|  | Penultimate | Penultimate | Last | Last | Next | Next
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
Sagittarius | $10^{9}$ | $2.7\times 10^{-1}$ | $-1.01$ | $4.9\times 10^{-8}$ | $-0.35$ | $1.3\times 10^{-1}$ | $0.03$
Hercules | $7.1\times 10^{6}$ | $8.4\times 10^{-8}$ | $-3.78$ | $2.4\times 10^{-3}$ | $-0.5$ | $2.4\times 10^{-3}$ | $3.18$
Leo II | $8.2\times 10^{6}$ | – | $-3.86$ | $1.6\times 10^{-3}$ | $-1.78$ | $3.2\times 10^{-3}$ | $2.31$
Segue 2 | $5.5\times 10^{5}$ | $6.2\times 10^{-4}$ | $-0.84$ | $8.5\times 10^{-4}$ | $-0.25$ | $6.3\times 10^{-5}$ | $0.28$
LMC | $1.4\times 10^{11}$ | $5.1\times 10^{-2}$ | $-7.63$ | – | $-2.67$ | $2.3\times 10^{-2}$ | $0.11$
SMC | $6.5\times 10^{9}$ | $2.8\times 10^{-5}$ | $-3.32$ | – | $-1.44$ | $8.7\times 10^{-6}$ | $0.22$
Draco I | $2.2\times 10^{7}$ | – | $-2.46$ | $9.9\times 10^{-5}$ | $-1.24$ | $7.1\times 10^{-6}$ | $0.24$
Bootes I | $10^{7}$ | $1.8\times 10^{-7}$ | $-1.67$ | $3.7\times 10^{-5}$ | $-0.35$ | – | $0.88$
Willman I | $4\times 10^{5}$ | $1.6\times 10^{-8}$ | $-0.66$ | $1.2\times 10^{-6}$ | $-0.21$ | $9.3\times 10^{-6}$ | $0.41$
Ursa Minor | $2\times 10^{7}$ | – | $-2.28$ | $1.7\times 10^{-5}$ | $-1.17$ | $2.6\times 10^{-6}$ | $0.29$
Ursa Major II | $4.9\times 10^{6}$ | $5.8\times 10^{-6}$ | $-2.12$ | $2.5\times 10^{-6}$ | $-0.09$ | – | $0.97$
Coma Berenices I | $1.2\times 10^{6}$ | $9.2\times 10^{-7}$ | $-2.58$ | $3.7\times 10^{-8}$ | $-0.25$ | – | $0.71$
Sculptor | $3.1\times 10^{7}$ | – | $-2.74$ | $3.4\times 10^{-8}$ | $-0.46$ | – | $1.48$
Table 1: Steady state response of the MW disk to encounters with satellites in
the collisionless limit, for the $(n,l,m)=(1,0,0)$ mode and for stars with
$I_{z}=I_{z,\odot}=h_{z}\sigma_{z,\odot}$ in the Solar neighborhood. We have
marginalized the response over $I_{R}$. Columns (1) and (2) list the name and
dynamical mass of each satellite. The latter is taken from the literature
(Simon & Geha, 2007; Bekki & Stanimirović, 2009; Łokas, 2009; Erkal et al.,
2019; Vasiliev & Belokurov, 2020), except for Sagittarius for which we adopt a
mass of $10^{9}\>{\rm M_{\odot}}$. Note that there is a discrepancy between
its estimated mass of $\sim 4\times 10^{8}\>{\rm M_{\odot}}$ (Vasiliev &
Belokurov, 2020) and the mass required ($10^{9}-10^{10}\>{\rm M_{\odot}}$) to
produce detectable phase-spiral signatures in N-body simulations (see for
example Bennett et al., 2022). Columns (3) and (4) respectively denote the
bending mode response assuming our fiducial MW parameters and the penultimate
disk-crossing time. Columns (5) and (6) indicate the same for the last disk-
crossing, while columns (7) and (8) show it for the next one. Only satellites
that induce a bending mode response, $f_{1,n=1}/f_{0}\geq 10^{-8}$, in at
least one of the three cases are shown. Any response weaker than $10^{-8}$ is
considered negligible and is indicated with a horizontal dash.

###  C.2 Special case: disk response for face-on impulsive encounters

The disk response in the general case, expressed by equation (38), depends on
several encounter parameters: $R_{\rm d}$, $\theta_{\rm P}$, $\phi_{\rm P}$,
and is complicated to evaluate. Therefore, as a sanity check, here we compute
the response as well as corresponding energy change for the special case of a
satellite undergoing an impulsive, perpendicular passage through the center of
the disk.

As shown in van den Bosch et al. (2018) (see also Banik & van den Bosch,
2021b), the total energy change due to a head-on encounter of velocity $v_{\rm
P}$ with a Plummer sphere of mass $M_{\rm P}$ and size $\varepsilon$ is given
by:

$\displaystyle\Delta E=4\pi\left({GM_{\rm p}\over v_{\rm
P}}\right)^{2}\int_{0}^{\infty}I^{2}_{0}(R)\Sigma(R){{\rm d}R\over R}$ (C11)

where

$\displaystyle I_{0}(R)=\int_{1}^{\infty}{M_{\rm P}(\zeta R)\over M_{\rm
p}}\,{{\rm d}\zeta\over\zeta^{2}(\zeta^{2}-1)^{1/2}}$ (C12)

Using that the enclosed mass profile of a Plummer sphere is given by $M_{\rm
P}(R)=M_{\rm P}R^{3}(R^{2}+\varepsilon^{2})^{-3/2}$, we have that
$I_{0}(R)=R^{2}/(R^{2}+\varepsilon^{2})$, which yields

$\displaystyle\Delta E=4\pi\left({GM_{\rm p}\over v_{\rm
P}}\right)^{2}\int_{0}^{\infty}\Sigma(R){R^{3}{\rm
d}R\over(R^{2}+\varepsilon^{2})^{2}}.$ (C13)

Now we compute the disk response to the face-on satellite encounter using
equations (20) and (C7-C10). For a perpendicular face-on impact through the
center of the disk we have $R_{\rm d}=0$ and $\theta_{\rm P}=0$, implying that
${\cal R}_{\rm c}$ becomes $R_{\rm c}$. The corresponding response is greatly
simplified. In the large time and small $I_{R}$ limit, it is given by equation
(20) with

$\displaystyle{\cal I}_{nlm}({\bf I},t)$ $\displaystyle\approx-\frac{2GM_{\rm
P}}{v_{\rm P}}\exp{\left[-i\Omega
t\right]}\,\delta_{m,0}\times\frac{1}{2\pi}\int_{0}^{2\pi}dw_{z}\exp{\left[-inw_{z}\right]}\exp{\left[i\frac{\Omega
z}{v_{\rm P}}\right]}$ $\displaystyle\times\frac{1}{2\pi}\int_{0}^{2\pi}{\rm
d}w_{R}\,\exp{\left[-ilw_{R}\right]}\,K_{0}\left[\frac{\left|\Omega\right|}{v_{\rm
P}}\sqrt{\varepsilon^{2}+{\left(R_{\rm
c}+\sqrt{\frac{2I_{R}}{\kappa}}\sin{w_{R}}\right)}^{2}}\right],$ (C14)

where the $\phi$ integral only leaves contribution from the axisymmetric $m=0$
mode. The $w_{R}$ integrand can be expanded as a Taylor series and the $w_{R}$
integral can be performed to yield the following leading order expression for
${\cal I}_{nlm}$:

$\displaystyle{\cal I}_{nlm}({\bf I},t)$ $\displaystyle\approx i\frac{GM_{\rm
P}}{v_{\rm P}}\exp{\left[-i\Omega
t\right]}\,\delta_{m,0}\,\left(\delta_{l,1}-\delta_{l,-1}\right)\times\frac{1}{2\pi}\int_{0}^{2\pi}dw_{z}\exp{\left[-inw_{z}\right]}\exp{\left[i\frac{\Omega
z}{v_{\rm P}}\right]}$
$\displaystyle\times\sqrt{\frac{2I_{R}}{\kappa}}\frac{R_{\rm
c}}{\sqrt{\varepsilon^{2}+R^{2}_{c}}}\,\frac{\left|\Omega\right|}{v_{\rm
P}}K^{\prime}_{0}\left[\frac{\left|\Omega\right|}{v_{\rm
P}}\sqrt{\varepsilon^{2}+R^{2}_{c}}\right].$ (C15)

In the impulsive limit, $v_{\rm P}\to\infty$, this becomes

$\displaystyle{\cal I}_{nlm}({\bf I},t)$ $\displaystyle\approx
i\,\delta_{n,0}\delta_{m,0}\left(\delta_{l,1}-\delta_{l,-1}\right)\frac{GM_{\rm
P}}{v_{\rm P}}\sqrt{\frac{2I_{R}}{\kappa}}\frac{R_{\rm
c}}{\varepsilon^{2}+R^{2}_{c}}\exp{\left[-i\,l\kappa\,t\right]},$ (C16)