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

TL;DR

This paper presents a linear perturbative framework to analyze phase spirals in inhomogeneous galactic disks, revealing how transient and resonant perturbations, especially from satellite encounters like Sagittarius, generate observable phase-space structures.

Contribution

It introduces a novel formalism for modeling phase mixing in stellar disks influenced by external perturbations, accounting for inhomogeneity and damping effects, advancing understanding of galactic phase-spiral formation.

Findings

Transient perturbations with timescales comparable to vertical oscillations trigger phase spirals.

Collisional damping causes super-exponential decay of phase-spiral amplitude.

Distant satellite encounters can induce stronger bending modes in the disk.

Abstract

We develop a linear perturbative formalism to compute the response of an inhomogeneous stellar disk embedded in a non-responsive dark matter halo to perturbations like bars, spiral arms and satellite galaxy encounters. 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 1-armed (2-armed) for vertically anti-symmetric (symmetric) bending (breathing) modes. Only transient perturbations with timescales ($\tau_{\mathrm{P}}$) comparable to the vertical oscillation period ($\tau_z \sim 1/\Omega_z$) trigger $z-v_z$ phase-spirals. Each $(n,l,m)$ mode of the response to impulsive ($\tau_{\mathrm{P}}<\tau=1/(n\Omega_z+l\Omega_r+m\Omega_\phi)$) perturbations is power law ($\sim \tau_{\mathrm{P}}/\tau$) suppressed, but that to adiabatic ($\tau_{\mathrm{P}}>\tau$) perturbations is exponentially weak ($\sim \exp{\left[-\left(\tau_{\mathrm{P}}/\tau\right)^\alpha\right]}$) except resonant ($\tau\to \infty$) modes. Slower ($\tau_{\mathrm{P}}>\tau_z$) perturbations, e.g., distant encounters with satellite galaxies, induce stronger bending modes. If the Gaia phase-spiral was triggered by a satellite, Sagittarius is the leading contender as it dominates the Solar neighborhood response of the Milky Way disk to satellite encounters. However, survival against collisional damping necessitates that the impact occurred within $\sim 0.6-0.7$ Gyr ago. We discuss how the detailed galactic potential dictates the phase-spiral shape: phase mixing occurs slower and phase-spirals are less wound in the outer disk and in presence of an ambient halo.

Figures

Tables

Table 1. 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 ) 𝑛 𝑙 𝑚 1 0 0 (n,l,m)=(1,0,0) mode and for stars with I z = I z , ⊙ = h z ​ σ z , ⊙ subscript 𝐼 𝑧 subscript 𝐼 𝑧 direct-product subscript ℎ 𝑧 subscript 𝜎 𝑧 direct-product I_{z}=I_{z,\odot}=h_{z}\sigma_{z,\odot} in the Solar neighborhood. We have marginalized the response over I R subscript 𝐼 𝑅 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 ​ M ⊙ superscript 10 9 subscript M direct-product 10^{9}\>{\rm M_{\odot}} . Note that there is a discrepancy between its estimated mass of ∼ 4 × 10 8 ​ M ⊙ similar-to absent 4 superscript 10 8 subscript M direct-product \sim 4\times 10^{8}\>{\rm M_{\odot}} (Vasiliev & Belokurov, 2020 ) and the mass required ( 10 9 − 10 10 ​ M ⊙ superscript 10 9 superscript 10 10 subscript M direct-product 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 ≥ 10 − 8 subscript 𝑓 1 𝑛 1 subscript 𝑓 0 superscript 10 8 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 superscript 10 8 10^{-8} is considered negligible and is indicated with a horizontal dash.

MW satellite	Mass	$f_{1,n=1}/f_0$	$t_{\mathrm{cross}}$	$f_{1,n=1}/f_0$	$t_{\mathrm{cross}}$	$f_{1,n=1}/f_0$	$t_{\mathrm{cross}}$
name	$\left(M_{\odot }\right)$		$\left(\mathrm{Gyr}\right)$		$\left(\mathrm{Gyr}\right)$		$\left(\mathrm{Gyr}\right)$
		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$

Equations

$$\frac{\partial f}{\partial t}+[f,H]=C[f],$$

$$C[f]=\frac{1}{2}\frac{\partial }{\partial {\xi }_i}\left(D_{ij}\frac{\partial f}{\partial {\xi }_j}\right),$$

$$[f_0,H_0]=C[f_0].$$

$$H=H_0+{\Phi }_{\mathrm{P}}(t)+{\Phi }_1(t),$$

$${\nabla }^2{\Phi }_1=4\pi G{\rho }_1.$$

$$f=f_0+f_1,$$

$$\frac{\partial f_1}{\partial t}+[f_1,H_0]+[f_0,{\Phi }_{\mathrm{P}}]+[f_0,{\Phi }_1]=C[f_1].$$

$${\Omega }_R=\frac{\partial H_0}{\partial I_R},{\Omega }_{\varphi }=\frac{\partial H_0}{\partial I_{\varphi }},{\Omega }_z=\frac{\partial H_0}{\partial I_z}.$$

$$\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 }_{\varphi }\frac{\partial f_1}{\partial w_{\varphi }}-\frac{\partial {\Phi }_{\mathrm{P}}}{\partial w_z}\frac{\partial f_0}{\partial I_z}-\frac{\partial {\Phi }_{\mathrm{P}}}{\partial w_R}\frac{\partial f_0}{\partial I_R}-\frac{\partial {\Phi }_{\mathrm{P}}}{\partial w_{\varphi }}\frac{\partial f_0}{\partial I_{\varphi }}$$

$$+\frac{D_z}{4I_z}\frac{{\partial }^2{f}_1}{\partial w_z^2}+\frac{D_R}{4I_R}\frac{{\partial }^2{f}_1}{\partial w_R^2}$$

$${\Phi }_{\mathrm{P}}\left(\mathbf{w},\mathbf{I},t\right)$$

$$f_1\left(\mathbf{w},\mathbf{I},t\right)$$

$$\frac{\partial f_{1nlm}}{\partial t}+i(n{\Omega }_z+l{\Omega }_R+m{\Omega }_{\varphi })f_{1nlm}$$

$$+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}.$$

$$f_{1nlm}(\mathbf{I},t)$$

$${\mathcal{I}}_{nlm}(\mathbf{I},t)$$

$${\mathcal{G}}_{nlm}(\mathbf{I},t-\tau )$$

$$\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],$$

$$f_0\approx \frac{\sqrt{2}}{{\pi }^{3/2}{\sigma }_z{h}_z}{\left(\frac{{\Omega }_{\varphi }\Sigma }{\kappa {\sigma }_R^2}\right)}_{R_{\mathrm{c}}}\exp \left[-\frac{\kappa I_R}{{\sigma }_R^2}\right]\exp \left[-\frac{E_z(I_z)}{{\sigma }_z^2}\right]\Theta (L_z),$$

$${\Phi }_0(R,z)={\Phi }_{\mathrm{d}}(R,z)+{\Phi }_{\mathrm{h}}(R,z).$$

$${\Phi }_z(R_{\mathrm{c}},z)={\Phi }_0(R_{\mathrm{c}},z)-{\Phi }_0(R_{\mathrm{c}},0).$$

$$I_z=\frac{1}{2\pi }\oint v_{z}dz=\frac{2}{\pi }{\int }_0^{z_{\max }}\sqrt{2[E_z-{\Phi }_z(R_{\mathrm{c}},z)]}dz,$$

$$T_z(R_{\mathrm{c}},I_z)=\oint \frac{dz}{v_z}=4{\int }_0^{z_{\max }}\frac{dz}{\sqrt{2\left[E_z(R_{\mathrm{c}},I_z)-{\Phi }_z(R_{\mathrm{c}},z)\right]}},$$

$$f_{1nlm}(\mathbf{I},t)$$

$$\times \left[\left\{\left(\frac{n{\Omega }_z}{{\sigma }_z^2}+\frac{l\kappa }{{\sigma }_R^2}\right)\left(\frac{{\Omega }_{\varphi }\Sigma }{\kappa }\right)-m\frac{\mathrm{d}}{\mathrm{d}{L}_z}\left(\frac{{\Omega }_{\varphi }\Sigma }{\kappa }\right)\right\}\Theta (L_z)-m\frac{{\Omega }_{\varphi }\Sigma }{\kappa }\delta (L_z)\right]{\mathcal{I}}_{nlm}(\mathbf{I},t).$$

$${\Phi }_{\mathrm{P}}(R,\varphi ,z)$$

$${\mathcal{M}}_j(t)$$

$${\mathcal{F}}_{\mathrm{o}}(z)$$

$${\mathcal{F}}_{\mathrm{e}}(z)$$

$${\mathcal{I}}_{nlm}(\mathbf{I},t)$$

$${\mathcal{P}}_{nlm}^{(j)}(\mathbf{I},t)$$

$$\times {\int }_0^{\infty }\mathrm{d}\tau \exp \left[-i{\Omega }_{\mathrm{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]{\mathcal{M}}_j(t-\tau ),$$

$${\Omega }_{\mathrm{res}}=n{\Omega }_z+l\kappa +m({\Omega }_{\varphi }-{\Omega }_{\mathrm{P}}).$$

$${\mathcal{P}}_{nlm}^{(j)}(\mathbf{I},t)$$

Full text

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:

[TABLE]

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):

[TABLE]

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),

[TABLE]

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

[TABLE]

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,

[TABLE]

The perturbed DF can be written as

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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}$:

[TABLE]

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}$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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,

[TABLE]

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}$):

[TABLE]

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:

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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:

[TABLE]

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,

[TABLE]

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

[TABLE]

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),

[TABLE]

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.,

[TABLE]

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

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}$.

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)

[TABLE]

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:

[TABLE]

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

[TABLE]

while that of a star is given by

[TABLE]

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

[TABLE]

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

[TABLE]

where $\Omega$ is given by

[TABLE]

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.

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:

[TABLE]

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

[TABLE]

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).

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

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,

[TABLE]

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:

[TABLE]

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:

[TABLE]

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.

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.

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

[TABLE]

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

Here $K_{0i}$ is defined as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

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.

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:

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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}$:

[TABLE]

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

[TABLE]

which can be substituted in equation (20) to yield

[TABLE]

with $f_{0}$ given by equation (15). Hence, the response is given by

[TABLE]

which shows that the satellite passage introduces a relative overdensity, $f_{1}\left({\bf w},{\bf I},t\right)/f_{0}({\bf I})$, that scales as $\sim R_{\rm c}/\left(\varepsilon^{2}+R^{2}_{c}\right)$, which increases from zero at the center, peaks at $R_{\rm c}=\varepsilon$, and asymptotes to zero again at large $R_{\rm c}$. The $\cos(w_{R}-\kappa t)$-term describes the radial epicyclic oscillations in the response.

To compute the energy change due to the impact, we note that ${\rm d}E/{\rm d}t=\partial E/\partial{\bf I}\cdot{\rm d}{\bf I}/{\rm d}t$, where $\partial E/\partial{\bf I}={\bf\Omega}=(\Omega_{z},\Omega_{R},\Omega_{\phi})$ and ${\rm d}{\bf I}/{\rm d}t=-\partial\Phi_{\rm P}/\partial{\bf w}$ from Hamilton’s equations of motion. Thus the total phase-averaged energy injected per unit phase-space can be obtained as follows:

[TABLE]

We can substitute the Fourier series expansions of $\Phi_{\rm P}$ and $f_{1}$ given in equations (10) in the above expression and integrate over ${\bf w}$ to obtain (Weinberg, 1994a, b)

[TABLE]

We can now substitute the form of $\Phi_{\rm P}$ for a Plummer perturber given in equation (37), with ${{\bf r}}_{\rm{{\bf P}}}$ and ${\bf r}$ given by equations (35) and (36). The time integral can thus be written as

[TABLE]

Using equations (B2) and (B4) to express $R$ and $z$ in terms of $({\bf w},{\bf I})$, and substituting the form for $f_{1nlm}({\bf I},t)$ from equation (C17), we can perform the above integrals over ${\bf w}$ and $t$. Substituting the result in equation (C20) we obtain

[TABLE]

The total energy, $\Delta E_{\rm tot}$, imparted into the disk by the impulsive satellite passage can be computed by integrating the above expression over ${\bf I}$ and ${\bf w}$ (which simply introduces a factor of ${\left(2\pi\right)}^{3}$ since $\left<\Delta E\left({\bf I}\right)\right>$ is already phase-averaged), using equation (15) and transforming from $L_{z}$ to $R_{\rm c}$ using the Jacobian ${\rm d}L_{z}/{\rm d}R_{\rm c}=R_{\rm c}\kappa^{2}/2\Omega_{\phi}$. This yields

[TABLE]

This is indeed the expression for $\Delta E_{\rm tot}$ derived under the impulse approximation given by equation (C13).