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

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 GMCs or DM substructure. The dynamics of stars in such a system is governed by the Fokker–Planck equation (FPE):

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+[f,H]=C[f],\end{eqnarray} \tag{ 1 }$$

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

$$\begin{eqnarray}&&C[f]=\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left({D}_{{ij}}\displaystyle \frac{\partial f}{\partial {\xi }_{j}}\right),\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{\xi }}=\left({\boldsymbol{q}},{\boldsymbol{p}}\right)$, with q and 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₀ and the corresponding DF be given by f₀, which satisfies the unperturbed FPE:

$$\begin{eqnarray}&&[{f}_{0},{H}_{0}]=C[{f}_{0}].\end{eqnarray} \tag{ 3 }$$

In the presence of a small time-dependent perturbation in the potential, Φ_P(t), the perturbed Hamiltonian can be written as

$$\begin{eqnarray}&&H={H}_{0}+{{\rm{\Phi }}}_{{\rm{P}}}(t)+{{\rm{\Phi }}}_{1}(t),\end{eqnarray} \tag{ 4 }$$

where Φ₁ is the gravitational potential related to the response density, ρ₁ = ∫f₁ d³ v , via the Poisson equation,

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{{\rm{\Phi }}}_{1}=4\pi G{\rho }_{1}.\end{eqnarray} \tag{ 5 }$$

The perturbed DF can be written as

$$\begin{eqnarray}&&f={f}_{0}+{f}_{1},\end{eqnarray} \tag{ 6 }$$

where f₁ is the linear-order perturbation in the DF. In the weak perturbation limit where linear perturbation theory holds, the time-evolution of f₁ is dictated by the following linearized form of the FPE:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{1}}{\partial t}+[{f}_{1},{H}_{0}]+[{f}_{0},{{\rm{\Phi }}}_{{\rm{P}}}]+[{f}_{0},{{\rm{\Phi }}}_{1}]=C[{f}_{1}].\end{eqnarray} \tag{ 7 }$$

Throughout this paper, we neglect the self-gravity of the disk response, which implies that we set the polarization term, [f₀, Φ₁] = 0. The implications of including self-gravity are discussed in Paper I.

The dynamics of a realistic disk galaxy like the MW is quasi-periodic, i.e., it can be characterized by oscillations in the azimuthal, radial, and vertical directions. In close proximity to the midplane and under radial epicyclic approximation, the Hamilton–Jacobi equation becomes separable, implying that all stars confined within a few vertical scale heights from the midplane of the disk are on regular, quasi-periodic orbits that are characterized by a radial action I_R , an azimuthal action I_ϕ , and a vertical action I_z . Hence, the motion of each star is characterized by three frequencies:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{R}=\frac{\partial {H}_{0}}{\partial {I}_{R}},\,\,\,{{\rm{\Omega }}}_{\phi }=\frac{\partial {H}_{0}}{\partial {I}_{\phi }},\,\,\,{{\rm{\Omega }}}_{z}=\frac{\partial {H}_{0}}{\partial {I}_{z}},\end{eqnarray} \tag{ 8 }$$

and three angles, ${w}_{i}(t)={w}_{i}(0)+{{\rm{\Omega }}}_{i}t$, with $i=R,\phi$ or $z$. This quasi-periodic nature of the orbits near the midplane is approximately preserved even in the presence of a (non-triaxial) DM halo, because this preserves the axisymmetry of the potential. Typically, as discussed in Section 2.3, 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

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1}}{\partial t}+{{\rm{\Omega }}}_{z}\displaystyle \frac{\partial {f}_{1}}{\partial {w}_{z}}+{{\rm{\Omega }}}_{R}\displaystyle \frac{\partial {f}_{1}}{\partial {w}_{R}}+{{\rm{\Omega }}}_{\phi }\displaystyle \frac{\partial {f}_{1}}{\partial {w}_{\phi }}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial {w}_{z}}\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{z}}\\ \quad -\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial {w}_{R}}\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{R}}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial {w}_{\phi }}\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{\phi }}={D}_{I}^{(z)}\displaystyle \frac{\partial }{\partial {I}_{z}}\left({I}_{z}\displaystyle \frac{\partial {f}_{1}}{\partial {I}_{z}}\right)\\ \quad +\displaystyle \frac{{D}_{I}^{(z)}}{4{I}_{z}}\displaystyle \frac{{\partial }^{2}{f}_{1}}{\partial {w}_{z}^{2}}+\displaystyle \frac{{D}_{I}^{(R)}}{4{I}_{R}}\displaystyle \frac{{\partial }^{2}{f}_{1}}{\partial {w}_{R}^{2}}.\end{array}\end{eqnarray} \tag{ 9 }$$

Here, we have performed several simplifications of the Fokker–Planck operator (see Appendix A for details). First, the diffusion coefficients are computed using the epicyclic approximation, i.e., small I_z and I_R , because only such stars are significantly affected by collisional scattering. In addition, following Tremaine et al. (2023), we consider ${D}_{I}^{(z)}$ and ${D}_{I}^{(R)}$ to be nearly constant, something that is implied by the age–velocity dispersion relation of the MW disk stars. Second, the I_R diffusion of the response f₁ is negligible 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. Third, following Binney & Lacey (1988), we have neglected diffusion in I_ϕ and w_ϕ , because the terms involving D_{ϕ ϕ} , D_{r ϕ} , and D_{ϕ z} are smaller than the I_z and I_R diffusion terms by factors of at least σ_R /v_c or σ_z /v_c , which are typically much smaller than unity (σ_R and σ_z are radial and vertical velocity dispersions, respectively, and v_c is the circular velocity along ϕ). 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 timescales much longer than that in actions—and hence it is comparatively less important.

Because the stars move along quasi-periodic orbits characterized by actions and angles, we can expand the perturbations, Φ_P and f₁, as discrete Fourier series in the angles as follows:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{\rm{P}}}\left({\boldsymbol{w}},{\boldsymbol{I}},t\right)\\ \,=\,\displaystyle \sum _{n=-\infty }^{\infty }\displaystyle \sum _{l=-\infty }^{\infty }\displaystyle \sum _{m=-\infty }^{\infty }\exp \left[i({{nw}}_{z}+{{lw}}_{R}+{{mw}}_{\phi })\right]\,{{\rm{\Phi }}}_{{nlm}}\left({\boldsymbol{I}},t\right),\\ {f}_{1}\left({\boldsymbol{w}},{\boldsymbol{I}},t\right)\\ \,=\,\displaystyle \sum _{n=-\infty }^{\infty }\displaystyle \sum _{l=-\infty }^{\infty }\displaystyle \sum _{m=-\infty }^{\infty }\exp \left[i({{nw}}_{z}+{{lw}}_{R}+{{mw}}_{\phi })\right]\,{f}_{1{nlm}}({\boldsymbol{I}},t),\end{array}\end{eqnarray} \tag{ 10 }$$

where w = (w_z , w_R , w_ϕ ) and I = (I_z , I_R , I_ϕ ). Substituting these Fourier expansions in Equation (9) yields the following differential equation for the evolution of f_1nlm :

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{1{nlm}}}{\partial t}+i(n{{\rm{\Omega }}}_{z}+l{{\rm{\Omega }}}_{R}+m{{\rm{\Omega }}}_{\phi }){f}_{1{nlm}}\\ \quad =i\left(n\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{z}}+l\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{R}}+m\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{\phi }}\right){{\rm{\Phi }}}_{{nlm}}\\ \qquad +{D}_{I}^{(z)}\displaystyle \frac{\partial }{\partial {I}_{z}}\left({I}_{z}\displaystyle \frac{\partial {f}_{1{nlm}}}{\partial {I}_{z}}\right)-\left[\displaystyle \frac{{n}^{2}{D}_{I}^{(z)}}{4{I}_{z}}+\displaystyle \frac{{l}^{2}{D}_{I}^{(R)}}{4{I}_{R}}\right]{f}_{1{nlm}}.\end{array}\end{eqnarray} \tag{ 11 }$$

This can be solved using the Green's function technique, with the initial condition, f_1nlm (t_i) = 0, to yield the following closed integral form for f_1nlm :

$$\begin{eqnarray}&&{f}_{1{nlm}}({\boldsymbol{I}},t)=i\left(n\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{z}}+l\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{R}}+m\displaystyle \frac{\partial {f}_{0}}{\partial {I}_{\phi }}\right){{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t).\end{eqnarray} \tag{ 12 }$$

Here, for ${D}_{I}^{(z)}\ll {\sigma }_{z}^{2}$ (σ_z is the vertical velocity dispersion), which is typically the case, ${{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t)$ can be approximately expressed as

$$\begin{eqnarray}&&{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t)\approx {\int }_{{t}_{{\rm{i}}}}^{t}d\tau \,{{ \mathcal G }}_{{nlm}}({\boldsymbol{I}},t-\tau ){{\rm{\Phi }}}_{{nlm}}({\boldsymbol{I}},\tau ).\end{eqnarray} \tag{ 13 }$$

Here, ${{ \mathcal G }}_{{nlm}}(t-\tau )$ is the Green's function (see Appendix A for detailed derivation), given by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal G }}_{{nlm}}({\boldsymbol{I}},t-\tau ) & \approx & \exp \left[-i(n{{\rm{\Omega }}}_{z}+l{{\rm{\Omega }}}_{R}+m{{\rm{\Omega }}}_{\phi })(t-\tau )\right]\\ & & \times \exp \left[-\left(\displaystyle \frac{{n}^{2}{D}_{I}^{(z)}}{4{I}_{z}}+\displaystyle \frac{{l}^{2}{D}_{I}^{(R)}}{4{I}_{R}}\right)\left(t-\tau \right)\right]\\ & & \times \,\exp \left[-\displaystyle \frac{{\left(n{{\rm{\Omega }}}_{z1}\right)}^{2}{D}_{I}^{(z)}{I}_{z}}{3}{\left(t-\tau \right)}^{3}\right],\end{array}\end{eqnarray} \tag{ 14 }$$

where Ω_z1 = ∂Ω_z /∂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Ω_z , lΩ_R , and mΩ_ϕ , which is being driven by an external time-dependent perturber potential, Φ_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 spatiotemporal behavior of the perturber potential, Φ_P, as well as the DF, f₀, of the unperturbed galaxy, which is addressed below.

Under the radial epicyclic approximation (small I_R ), the unperturbed DF, f₀, 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

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0}\approx \displaystyle \frac{\sqrt{2}}{{\pi }^{3/2}\,{\sigma }_{z}{h}_{z}}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{\phi }{\rm{\Sigma }}}{\kappa \,{\sigma }_{R}^{2}}\right)}_{{R}_{{\rm{c}}}}\\ \,\times \,\exp \left[-\displaystyle \frac{\kappa {I}_{R}}{{\sigma }_{R}^{2}}\right]\,\exp \left[-\displaystyle \frac{{E}_{z}({I}_{z})}{{\sigma }_{z}^{2}}\right]\,{\rm{\Theta }}({L}_{z}).\end{array}\end{eqnarray} \tag{ 15 }$$

The vertical structure of this disk is isothermal, while the radial profile is pseudo-isothermal. ³ Here, ${\rm{\Sigma }}={\rm{\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_ϕ , R_c = R_c(L_z ) is the guiding radius, Ω_ϕ is the circular frequency, and $\kappa =\kappa ({R}_{{\rm{c}}})={\mathrm{lim}}_{{I}_{R}\to 0}{{\rm{\Omega }}}_{R}$ is the radial epicyclic frequency (Binney & Tremaine 1987). Θ(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, ρ(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 (sech²) profile with thickness h_z (Equation (B3)). As shown by Smith et al. (2015), this density profile is accurately approximated by a sum of three Miyamoto & Nagai (1975) disks, ⁴ which has a simple, analytical form for the associated potential. Therefore, we use this 3MN approximation for our disk throughout this work, because this drastically simplifies the computation of orbital frequencies. For the purpose of computing the disk response, we assume typical MW-like parameters for the various quantities, i.e., R_⊙ = 8 kpc, disk mass M_d = 5 × 10¹⁰ M_⊙, h_R = 2.2 kpc, σ_R (R_⊙) = σ_R,⊙ = 35 km s⁻¹, h_z = 0.4 kpc, and ${\sigma }_{z}({R}_{\odot })={\sigma }_{z,\odot }=\sqrt{2\pi {{Gh}}_{z}{\rm{\Sigma }}({R}_{\odot })}=23\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (McMillan 2011; Bovy & Rix 2013). For this set of parameters, the Toomre Q = σ_R κ/(3.36GΣ) for the stellar disk turns out to be 2.26 in the Solar neighborhood, indicating that the disk is gravitationally stable. The isothermal vertical distribution of disk stars adopted here ignores the potential influence of gas near the midplane, which may increase the shear, $\left|d{{\rm{\Omega }}}_{z}/d\mathrm{ln}{I}_{z}\right|$, for small I_z . This may affect the shape of the phase spiral, making it more tightly wound, but should not substantially impact its amplitude.

The disk is assumed to be embedded in an extended DM halo characterized by a spherical Navarro–Frenk–White (NFW; Navarro et al. 1997) density profile, with virial mass M_vir, concentration c, scale radius r_s , and the corresponding potential Φ_h given by Equation (B6). We adopt M_vir = 9.78 × 10¹¹ M_⊙, 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.,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{0}(R,z)={{\rm{\Phi }}}_{{\rm{d}}}(R,z)+{{\rm{\Phi }}}_{{\rm{h}}}(R,z).\end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{z}({R}_{{\rm{c}}},z)={{\rm{\Phi }}}_{0}({R}_{{\rm{c}}},z)-{{\rm{\Phi }}}_{0}({R}_{{\rm{c}}},0).\end{eqnarray} \tag{ 17 }$$

The vertical action, I_z , can be obtained from E_z as follows:

$$\begin{eqnarray}&&{I}_{z}=\displaystyle \frac{1}{2\pi }\oint {v}_{z}\,{dz}=\displaystyle \frac{2}{\pi }{\int }_{0}^{{z}_{\max }}\sqrt{2[{E}_{z}-{{\rm{\Phi }}}_{z}({R}_{{\rm{c}}},z)]}\,{dz},\end{eqnarray} \tag{ 18 }$$

where ${{\rm{\Phi }}}_{z}({R}_{{\rm{c}}},{z}_{\max })={E}_{z}$. This implicit equation can be inverted to obtain E_z (R_c, I_z ). The time period of vertical oscillation can then be obtained using

$$\begin{eqnarray}&&{T}_{z}({R}_{{\rm{c}}},{I}_{z})=\oint \displaystyle \frac{{dz}}{{v}_{z}}=4{\int }_{0}^{{z}_{\max }}\displaystyle \frac{{dz}}{\sqrt{2\left[{E}_{z}({R}_{{\rm{c}}},{I}_{z})-{{\rm{\Phi }}}_{z}({R}_{{\rm{c}}},z)\right]}},\end{eqnarray} \tag{ 19 }$$

which yields the vertical frequency, Ω_z (R_c, I_z ) = 2π/T_z (R_c, I_z ). The tightness of a vertical phase spiral at a given I_z increases for faster phase mixing, i.e., increasing value of the shear, $\left|d{{\rm{\Omega }}}_{z}/d\mathrm{ln}{I}_{z}\right|$, at that I_z . The phase-mixing timescale, τ_ϕ , is given by the inverse of the shear:

$$\begin{eqnarray}&&{\tau }_{\phi }=\left|\displaystyle \frac{d\mathrm{ln}{I}_{z}}{d{{\rm{\Omega }}}_{z}}\right|.\end{eqnarray} \tag{ 20 }$$

Figure 1 plots Ω_z and τ_ϕ as a function of I_z in the left and right panels, respectively, for four different cases: (i) MW disk plus halo at R_c = 8 kpc (solid red line), (ii) only the MW disk at R_c = 8 kpc (dashed red line), (iii) MW disk plus halo at R_c = 12 kpc (dotted–dashed blue line), and (iv) only the MW disk at R_c = 12 kpc (dotted blue line). The presence of the DM halo deepens and steepens the potential well and thus increases the vertical frequencies. In the presence of the halo and farther out in the disk, the phase-mixing timescale, τ_ϕ , is longer, i.e., phase mixing occurs more slowly, and the phase spiral is more loosely wound. The sensitivity of the tightness of the phase spiral to the detailed MW potential implies that the phase spiral can be used to constrain it.

Zoom In Zoom Out Reset image size

Substituting the expression of f₀ for a MW-like disk given by Equation (15) in Equation (12), we obtain the following integral form for f_1nlm :

$$\begin{eqnarray}\begin{array}{rcl}{f}_{1{nlm}}({\boldsymbol{I}},t) & \approx & -\displaystyle \frac{2i}{\pi {\sigma }_{R}^{2}}\,\displaystyle \frac{1}{\sqrt{2\pi }{h}_{z}{\sigma }_{z}}\exp \left[-\displaystyle \frac{\kappa {I}_{R}}{{\sigma }_{R}^{2}}\right]\exp \left[-\displaystyle \frac{{E}_{z}({I}_{z})}{{\sigma }_{z}^{2}}\right]\\ & & \times \left[\left\{\left(\displaystyle \frac{n{{\rm{\Omega }}}_{z}}{{\sigma }_{z}^{2}}+\displaystyle \frac{l\kappa }{{\sigma }_{R}^{2}}\right)\left(\displaystyle \frac{{{\rm{\Omega }}}_{\phi }{\rm{\Sigma }}}{\kappa }\right)-m\displaystyle \frac{d}{{{dL}}_{z}}\left(\displaystyle \frac{{{\rm{\Omega }}}_{\phi }{\rm{\Sigma }}}{\kappa }\right)\right\}\right.\\ & & \left.\times {\rm{\Theta }}({L}_{z})-m\displaystyle \frac{{{\rm{\Omega }}}_{\phi }{\rm{\Sigma }}}{\kappa }\delta ({L}_{z})\right]{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t).\end{array}\end{eqnarray} \tag{ 21 }$$

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. Second, for the sake of simplicity, we consider the ambient DM halo to be nonresponsive, and therefore we 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 spatiotemporal nature of the perturbing potential dictates the disk response. In this paper, we explore two different types of perturbation to which realistic disk 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/bars. 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.

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 (C5) in Equation (21) and performing the τ integration with the initial time, t_i → − ∞ . This yields the modal response, f_1nlm (Equation (21)), with ${{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t)$ given by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & = & \alpha \,{{\rm{\Psi }}}_{{nlm}}^{({\rm{o}})}({\boldsymbol{I}}){{ \mathcal P }}_{{nlm}}^{({\rm{o}})}({\boldsymbol{I}},t)\\ & & +{{\rm{\Psi }}}_{{nlm}}^{({\rm{e}})}({\boldsymbol{I}}){{ \mathcal P }}_{{nlm}}^{({\rm{e}})}({\boldsymbol{I}},t),\end{array}\end{eqnarray} \tag{ 25 }$$

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

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal P }}_{{nlm}}^{(j)}({\boldsymbol{I}},t) & = & \exp \left[-i\,m{{\rm{\Omega }}}_{{\rm{P}}}\,t\right]\\ & & \times {\int }_{0}^{\infty }d\tau \exp \left[-i\,{{\rm{\Omega }}}_{\mathrm{res}}\,\tau \right]\\ & & \times \,\exp \left[-\left(\frac{{n}^{2}{D}_{I}^{(z)}}{4{I}_{z}}+\frac{{l}^{2}{D}_{I}^{(R)}}{4{I}_{R}}\right)\tau \right]\\ & & \times \,\exp \left[-\frac{{\left(n{{\rm{\Omega }}}_{z1}\right)}^{2}{D}_{I}^{(z)}{I}_{z}}{3}{\tau }^{3}\right]\,{{ \mathcal M }}_{j}(t-\tau ),\end{array}\end{eqnarray} \tag{ 26 }$$

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

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{res}}=n{{\rm{\Omega }}}_{z}+l\kappa +m({{\rm{\Omega }}}_{\phi }-{{\rm{\Omega }}}_{{\rm{P}}}).\end{eqnarray} \tag{ 27 }$$

3.1.1. Collisionless Limit

First, we examine the response in the limit of zero diffusion, i.e., ${D}_{I}^{(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 ${{ \mathcal M }}_{j}(t)$ is given by the first of Equations (23). In this case,

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal P }}_{{nlm}}^{(j)}({\boldsymbol{I}},t) & = & \displaystyle \frac{1}{2\,{\omega }_{j}}\,\exp \left[-\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{res}}^{2}}{4{\omega }_{j}^{2}}\right]\left[1+\mathrm{erf}\left({\omega }_{j}t-i\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{res}}}{2{\omega }_{j}}\right)\right]\\ & & \times \,\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi })t\right]\\ & \mathop{\to }\limits^{t\to \infty } & \displaystyle \frac{1}{{\omega }_{j}}\,\exp \left[-\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{res}}^{2}}{4{\omega }_{j}^{2}}\right]\\ & & \times \,\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi })t\right].\end{array}\end{eqnarray} \tag{ 28 }$$

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

The left-hand panel of Figure 2 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, ω_j (j = o and 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 Σ_P = 5.5 M_⊙ pc⁻², Ω_P = 12 km s⁻¹ kpc⁻¹, ${k}_{z}^{(o)}={k}_{z}^{(e)}=1\,\mathrm{kpc}$, k_R = 10 kpc, and I_z = I_z,⊙ ≡ h_z σ_z,⊙ = 9.2 km s⁻¹, and marginalize the response over I_R . We set α = 1, implying equal maximum strengths for the bending and breathing modes. As evident from this figure, and also from Equation (28), the long-term strength of the disk response (after the initial transients have died out like ${e}^{-{\omega }_{j}^{2}{t}^{2}}$) scales as ∼1/ω_j in the impulsive (large ω_j ) limit, but is superexponentially suppressed ($\sim \exp \left[-{{\rm{\Omega }}}_{\mathrm{res}}^{2}/4{\omega }_{j}^{2}\right]$) in the adiabatic (small ω_j ) limit away from resonances, i.e., for ${{\rm{\Omega }}}_{\mathrm{res}}\ne 0$. The adiabatic suppression scales differently with ω_j for other functional forms of ${{ \mathcal M }}_{j}(t)$; e.g., for ${{ \mathcal M }}_{j}(t)=1/\sqrt{1+{\omega }_{j}^{2}{t}^{2}}$, the response strength is exponentially suppressed $(\sim \exp [-{{\rm{\Omega }}}_{\mathrm{res}}/{\omega }_{j}])$. However, the response of resonant modes (${{\rm{\Omega }}}_{\mathrm{res}}=0$) does not undergo adiabatic suppression, and it scales as ∼1/ω_j throughout, becoming nonlinear in the adiabatic regime. Because 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 ω_j (Weinberg 1994a, 1994b).

Zoom In Zoom Out Reset image size

The sinusoidal factor, $\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi })t\right]$, in ${{ \mathcal P }}_{{nlm}}^{(j)}$ describes the oscillations of stars with three different frequencies, Ω_z , κ, and Ω_ϕ , along the vertical, radial, and azimuthal directions, respectively. Due to the dependence of these frequencies on the actions, that of Ω_z on I_z and of κ and Ω_ϕ on I_ϕ = 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 \mbox{--}\dot{\phi }$ phase spaces, respectively. As is evident from Equation (26), ${{ \mathcal P }}_{{nlm}}^{(j)}\sim \exp \left[-{im}{{\rm{\Omega }}}_{{\rm{P}}}t\right]$ in the adiabatic limit (ω_j → 0); hence, in this limit, the sinusoidal factor $\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\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 because Ω_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. Because Ω_z , Ω_ϕ , and Ω_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 ΔL_z around L_z , damps out as follows:

$$\begin{eqnarray}\begin{array}{rcl}\left\langle {f}_{1{nlm}}\right\rangle ({\boldsymbol{I}},t) & = & \displaystyle \frac{1}{{\rm{\Delta }}{L}_{z}}{\int }_{{L}_{z}-{\rm{\Delta }}{L}_{z}/2}^{{L}_{z}+{\rm{\Delta }}{L}_{z}/2}{{dL}}_{z}\,{f}_{1{nlm}}({\boldsymbol{I}},t)\\ & \approx & \displaystyle \frac{\sin \left[\left(\tfrac{\partial }{\partial {L}_{z}}\left(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi }\right)\right)\tfrac{{\rm{\Delta }}{L}_{z}}{2}t\right]}{\left(\tfrac{\partial }{\partial {L}_{z}}\left(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi }\right)\right)\tfrac{{\rm{\Delta }}{L}_{z}}{2}t}\,{f}_{1{nlm}}({\boldsymbol{I}},t).\end{array}\end{eqnarray} \tag{ 29 }$$

Because the frequencies vary with L_z , marginalizing over L_z mixes phase spirals that differ slightly in phases, giving way to a ∼1/t damping accompanied by a beat-like modulation with a characteristic lateral mixing timescale:

$$\begin{eqnarray}&&{\tau }_{{\rm{D}}}^{(\mathrm{LM})}=\displaystyle \frac{1}{\left(\tfrac{\partial }{\partial {L}_{z}}\left(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi }\right)\right)\tfrac{{\rm{\Delta }}{L}_{z}}{2}}.\end{eqnarray} \tag{ 30 }$$

This explains why the density contrast of the Gaia phase spiral is enhanced upon color coding by v_ϕ , 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 ∼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

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{res}}=n{{\rm{\Omega }}}_{z}+l\kappa +m({{\rm{\Omega }}}_{\phi }-{{\rm{\Omega }}}_{{\rm{P}}})=0,\end{eqnarray} \tag{ 31 }$$

do not phase mix away. Hence, parts of the phase space closer to a resonance take longer to phase mix away. Moreover, as is manifest from the adiabatic suppression factor, $\exp [-{{\rm{\Omega }}}_{\mathrm{res}}^{2}/4{\omega }_{j}^{2}]$, the near-resonant modes with ${{\rm{\Omega }}}_{\mathrm{res}}\ll 2{\omega }_{j}$ have much larger amplitude than those with ${{\rm{\Omega }}}_{\mathrm{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 corotation resonance (n, l, m) = (0, 0, m), the Lindblad resonances (0, ±1, ±2), the ultraharmonic resonances (0, ±1, ±4), and so on. For thin disks with h_z ≪ h_R , the vertical degrees of freedom are generally not in resonance with the radial or azimuthal ones, because Ω_z is much larger than Ω_ϕ or κ. Hence, the vertical oscillation modes (n ≠ 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 (Ω_z = 2Ω_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, ω_j , and the resonant frequency, ${{\rm{\Omega }}}_{\mathrm{res}}$. The right panel of Figure 2 shows the breathing-to-bending ratio, f_1,200/f_1,100, as a function of ω_e and ω_o, with blue (yellow) shades indicating low (high) values. In general, the breathing-to-bending ratio rises steeply and falls gradually with ω_e at fixed ω_o, while the trend is reversed as a function of ω_o at fixed ω_e, resulting in a saddle point at (ω_e, ω_o) ≈ (9, 7). This owes to the superexponential suppression in the adiabatic (${\omega }_{j}\ll {{\rm{\Omega }}}_{\mathrm{res}}$) limit and the power-law suppression in the impulsive (${\omega }_{j}\gg {{\rm{\Omega }}}_{\mathrm{res}}$) limit. Along the ω_o = ω_e line, the bending (breathing) modes dominate in the adiabatic (impulsive) limit, which is also evident from the left panel of Figure 2. All this suggests that bending modes dominate over breathing modes when (i) the antisymmetric perturbation is more impulsive, i.e., evolves faster than the symmetric one, or (ii) both symmetric and antisymmetric perturbations occur over comparable timescales but more slowly 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 ${{ \mathcal M }}_{j}(t)$ is given by the second form of Equation (23). In this case, as shown by Equation (19) of Banik & van den Bosch (2021a):

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal P }}_{{nlm}}^{(j)}({\boldsymbol{I}},t)=\displaystyle \frac{\exp \left[{\gamma }_{j}t\right]\exp \left[-{im}{{\rm{\Omega }}}_{{\rm{P}}}t\right]}{{\gamma }_{j}+i{{\rm{\Omega }}}_{\mathrm{res}}}\left[1-{\rm{\Theta }}(t)\right]\\ \,+\,i\left[\displaystyle \frac{{\gamma }_{j}\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi })t\right]}{{{\rm{\Omega }}}_{\mathrm{res}}({\gamma }_{j}+i{{\rm{\Omega }}}_{\mathrm{res}})}-\displaystyle \frac{\exp \left[-{im}{{\rm{\Omega }}}_{{\rm{P}}}t\right]}{{{\rm{\Omega }}}_{\mathrm{res}}}\right]{\rm{\Theta }}(t).\end{array}\end{eqnarray} \tag{ 32 }$$

Up to t = 0, when the perturber amplitude stops growing, the response from all modes oscillates with the pattern speed Ω_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 Ω_P.

For a slowly growing spiral/bar, i.e., in the "adiabatic growth" limit (γ → 0), the entire disk oscillates at the driving frequency, Ω_P, i.e.,

$$\begin{eqnarray}&&{{ \mathcal P }}_{{nlm}}^{(j)}({\boldsymbol{I}},t)\mathop{\to }\limits^{{\gamma }_{j}\to 0}\exp \left[-{im}{{\rm{\Omega }}}_{{\rm{P}}}t\right]\left[\pi \delta ({{\rm{\Omega }}}_{\mathrm{res}})-\displaystyle \frac{i}{{{\rm{\Omega }}}_{\mathrm{res}}}\right].\end{eqnarray} \tag{ 33 }$$

This has two major implications. First of all, because all stars, both resonant and nonresonant, are driven at the pattern speed of the perturbing spiral/bar, transient phase mixing does not occur and thus no phase spiral arises. Second, the response is dominated by the resonances, ${{\rm{\Omega }}}_{\mathrm{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; Banik & van den Bosch 2022; Chiba & Schönrich 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 corotation 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 are 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}_{I}^{(z)}$ and ${D}_{I}^{(R)}$. Following Tremaine et al. (2023), we assume that the disk stars have gained their mean vertical and radial actions over the age of the disk, T_disk = 10 Gyr, due to collisional heating, which implies that ${D}_{a}=\left\langle {I}_{a}\right\rangle /{T}_{\mathrm{disk}}$, where a is either z or R and $\left\langle {I}_{a}\right\rangle =\int {{dI}}_{a}\,{I}_{a}\,{f}_{0}\,/\int {{dI}}_{a}\,{f}_{0}$.

For a transient bar/spiral with pulse frequency ω_j , ${{ \mathcal P }}_{{nlm}}^{(j)}$ is given by Equation (26). In the impulsive limit (ω_j → ∞ ), we have that ${{ \mathcal M }}_{j}(t-\tau )\to {\omega }_{j}\delta (t-\tau )$. Upon absorbing ω_j in the prefactor, the expression for ${{ \mathcal P }}_{{nlm}}^{(j)}$ then simplifies to

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal P }}_{{nlm}}^{(j)}({\boldsymbol{I}},t)\approx {\rm{\Theta }}(t)\exp \left[-i(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi })t\right]\\ \,\times \,\exp \left[-\left(\displaystyle \frac{{n}^{2}{D}_{I}^{(z)}}{4{I}_{z}}+\displaystyle \frac{{l}^{2}{D}_{I}^{(R)}}{4{I}_{R}}\right)t\right]\,\exp \left[-\displaystyle \frac{{\left(n{{\rm{\Omega }}}_{z1}\right)}^{2}{D}_{I}^{(z)}{I}_{z}}{3}{t}^{3}\right].\end{array}\end{eqnarray} \tag{ 34 }$$

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:

$$\begin{eqnarray}&&{\tau }_{{\rm{D}}}^{({\rm{w}})}={\left[\displaystyle \frac{{n}^{2}{D}_{I}^{(z)}}{4{I}_{z}}+\displaystyle \frac{{l}^{2}{D}_{I}^{(R)}}{4{I}_{R}}\right]}^{-1},\quad {\tau }_{{\rm{D}}}^{({\rm{I}})}={\left[\displaystyle \frac{3}{{\left(n{{\rm{\Omega }}}_{z1}\right)}^{2}{D}_{I}^{(z)}{I}_{z}}\right]}^{1/3}.\end{eqnarray} \tag{ 35 }$$

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., Ω_z1 ≠ 0). The left (right) panel of Figure 3 plots the diffusion timescale, ${\tau }_{{\rm{D}}}^{({\rm{I}})}$, as a function of I_z (R_c) for three different values of R_c (I_z ) as indicated. We note that ${\tau }_{{\rm{D}}}^{({\rm{I}})}$ diverges in the small I_z limit, attains a minimum around I_z ∼ 0.2 − 0.5 h_z σ_z,⊙, and increases as ${I}_{z}^{\beta }$ with β < 1 at large I_z . As a function of R_c, ${\tau }_{{\rm{D}}}^{({\rm{I}})}$ shows an approximately exponential rise. This is owed to the fact that $\left\langle {D}_{I}^{(z)}\right\rangle \sim {h}_{z}\,{\sigma }_{z}({R}_{{\rm{c}}})/{T}_{\mathrm{disk}}\,\sim \exp \left[-{R}_{{\rm{c}}}/2{h}_{R}\right]/{T}_{\mathrm{disk}}$ for the 3MN profile adopted for the MW disk. At R_c = R_⊙ = 8 kpc, ${\tau }_{{\rm{D}}}^{(I)}\sim 0.6\mbox{--}0.7\,\mathrm{Gyr}$, in agreement with Tremaine et al. (2023). Hence, we see that collisional diffusion in action space is fairly efficient, and thus that phase spirals are short-lived features.

Zoom In Zoom Out Reset image size

Figure 4 plots the amplitude of the disk response (for I_R = 0) to a transient spiral of pulse frequency, ω_o = ω_e = 0.5σ_z,⊙/h_z , computed using Equations (21), (25), and (26), 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_c, while the columns denote different values of ${I}_{z}/({h}_{z}{\sigma }_{z}({R}_{c}))$ 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 gray 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 (28) 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. We note that the collisional damping is faster for smaller R_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}})}$.

Zoom In Zoom Out Reset image size

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

It is instructive to study the two extreme cases of encounter speed, the impulsive limit (large v_P) and the adiabatic limit (small v_P). Using the asymptotic form of the K₀ Bessel function that appears in Equation (39), we obtain the following approximate asymptotic behavior of f_1nlm at large time:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{1{nlm}} & \sim & \displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[-i{\rm{\Omega }}t\right]\\ & & \times \left\{\begin{array}{ll}1, & {v}_{{\rm{P}}}\to \infty \\ \sqrt{{v}_{{\rm{P}}}/{\rm{\Omega }}b}\,\exp \left[-{\rm{\Omega }}b/{v}_{{\rm{P}}}\right], & {v}_{{\rm{P}}}\to 0,\end{array}\right.\end{array}\end{eqnarray} \tag{ 41 }$$

where b is the impact parameter of the encounter, defined as the perpendicular distance of the nearest star on the midplane from the satellite's (straight) orbit, and expressed as

$$\begin{eqnarray}&&b=\left|{R}_{{\rm{d}}}-{R}_{{\rm{c}}}\right|\,\sqrt{1-{\sin }^{2}{\theta }_{{\rm{P}}}{\cos }^{2}{\phi }_{{\rm{P}}}}.\end{eqnarray} \tag{ 42 }$$

It is clear from these limits that the disk response is most pronounced for intermediate velocities, v_P ∼ Ωb. For impulsive encounters, the response is suppressed as a power law in v_P, whereas in the adiabatic limit, the response is exponentially suppressed, except at resonances, Ω = nΩ_z + l κ + mΩ_ϕ = 0. In this limit, far from the resonances, the perturbation timescale, b/v_P, is much larger than Ω⁻¹, 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, 1994b; Gnedin & Ostriker 1999).

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 to satellite encounters that occurred in the past few hundred million years, which we may expect to have triggered phase spirals that 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 A2; 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 bulge, ⁵ using a second-order leapfrog integrator. For each individual orbit, we record the times, t_cross, and the galactocentric radii, R_d, corresponding to disk crossings. We also register the corresponding impact velocities, ${v}_{{\rm{P}}}=\sqrt{{v}_{z}^{2}+{v}_{R}^{2}+{v}_{\phi }^{2}}$, 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 (D7) and compute the disk response (integrated over I_R ) following the satellite encounter, using Equations (21) and (D7). Results are summarized in Table 1 of Appendix D.1. Figure 6 plots the amplitude of the solar neighborhood (for which R_c(L_z ) = R_⊙ = 8 kpc) bending mode response, f_1,n=1/f₀ (top panel), and breathing-to-bending ratio, f_1,n=2/f_1,n=1 (bottom panel), as a function of t_cross. Here, we only show the responses for (l, m) = (0, 0) modes and consider stars with I_z = I_z,⊙ = 9.2 kpc km s⁻¹.

Zoom In Zoom Out Reset image size

It is noteworthy that the responses in the realistic MW disk computed here are ∼1–2 orders of magnitude larger than those evaluated for the isothermal slab model shown in Figure 7 of Paper I. This is owed 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 Figure 6, 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}\gt {\delta }_{\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 insufficient 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₀, 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 ∼1 Gyr ago, triggered a strong response of f_1,n=1/f₀ ∼ 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 Myr ago, triggered a very weak, adiabatically suppressed response (∼5 × 10⁻⁸) that falls below the lower limit of Figure 6. Its next disk crossing in about 30 Myr is estimated to trigger a strong response with f_1,n=1/f₀ ∼ 0.1. Besides Sgr, the satellites that excite a detectable response, ${f}_{1,n=1}/{f}_{0}\gt {\delta }_{\min }$ are Hercules, Segue 2, Leo II, and the LMC. The imminent crossing of LMC is estimated to trigger f_1,n=1/f₀ ∼ 2 × 10⁻², which is an order of magnitude below Sgr. Only for I_z /I_z,⊙ ≳ 4.5 (${z}_{\max }\gtrsim 3.4{h}_{z}$) does the LMC response dominate 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; Bland-Hawthorn et al. 2019; Darling & Widrow 2019a; Bland-Hawthorn & Tepper-García 2021; Hunt et al. 2021).

We caution that the above estimates of the disk response are obtained under the assumption of nearly straight-line orbits of the satellites. This assumption is well-justified as long as the satellite orbit is sufficiently eccentric. However, realistic orbits with low eccentricities can trigger a substantially different disk response. Also, there is a large uncertainty in Sgr's orbit (as well as that of the other satellites), due to the uncertainty in its observed phase-space coordinates and the MW parameters. These can explain the apparent discrepancy between our results and those of Bennett et al. (2022): they find in their N-body simulations that Sgr triggers a weaker response in the Solar neighborhood than what we compute in this paper. This is probably attributable to the fact that Sgr's orbit induces an adiabatic perturbation in the solar neighborhood in their simulations. To integrate the satellite orbits for our analysis, we pick the median values of the observed phase-space coordinates quoted by Riley et al. (2019) and Vasiliev & Belokurov (2020) as the initial conditions. Other values within the margin of error may lead to substantially different orbits, implying different R_d, v_P, θ_P and ϕ_P, and therefore different responses. Moreover, some of these orbits can be mildly eccentric, for which the straight-line orbit approximation adopted in this paper is not well-justified. A perturbative analysis of the disk response to satellites along general orbits is beyond the scope of this paper and deserves future investigation. We also ignore the effect of dynamical friction on the satellite orbits and thus on the resulting disk response. 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 × 10¹² M_⊙ leaves the relative amplitudes of the satellite responses (in the collisionless limit) nearly unchanged, but it 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 ∼600 Myr ago (as opposed to 1 Gyr ago in the fiducial case), which is closer to the winding time of ∼500 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 superexponentially over a timescale that is ∼0.6–0.7 Gyr in the solar neighborhood (see Section 3.1.2). This would almost completely wash away the response to any satellite encounter that occurred ≳1 Gyr ago. For example, the present-day response to the last pericentric passage of Leo II that occurred ∼1.8 Gyr ago would be completely erased. If the last disk crossing of Sgr that induced a strong response occurred ∼1 Gyr ago, as in the fiducial MW model, the response would have been damped out by ∼2 orders of magnitude by today, meaning Sgr is 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 × 10¹² M_⊙ implies a Sgr crossing time of ∼0.6 Gyr instead of 1 Gyr. In this case, the response would only have been damped by a factor of ∼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 ∼0.6–0.7 Gyr from the present day.

Having computed the MW disk response to its satellites, we now investigate the sensitivity of the response to the various encounter parameters. In Figure 7, we plot the amplitude of the solar neighborhood response, f_1,nlm /f₀ (marginalized over I_R ), as a function of the impact velocity, v_P (in units of the circular velocity at R_c = R_⊙), 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 , θ_P, and ϕ_P, respectively, assuming the fiducial parameters to be those for Sgr (mass M_P = 10⁹ M_⊙, scale radius ε = 1.6 kpc) during its penultimate disk crossing (most relevant for the Gaia snail), i.e., impact radius R_d = 17 kpc, impact velocity v_P = 340 km s⁻¹, and angles of impact, θ_P = 21° and ϕ_P = 150°. In Figure 8, we plot the bending and breathing mode response amplitudes (in the solar neighborhood) as a function of v_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. Figure 9 shows the ratio of the bending and breathing response amplitudes as a function of v_P for the dominant mode (l, m) = (0, − 2). Different lines indicate breathing-to-bending ratios for different values of θ_P, while the left and right columns, respectively, indicate the ratios observed at R_c = 8 and 12 kpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 7 and 8, it is evident that, as shown in Equation (41), the disk response is suppressed like a power law ($\sim {v}_{{\rm{P}}}^{-1}$) in the high-velocity/impulsive limit and exponentially ($\sim \exp \left[-{\rm{\Omega }}b/{v}_{{\rm{P}}}\right]$) suppressed in the low-velocity/adiabatic limit. The response is the strongest for intermediate velocities, v_P ∼ 2−3 v_circ(R_⊙), where the frequencies of the vertical, radial, and azimuthal oscillations of the stars are nearly commensurate with the encounter frequency, ${v}_{{\rm{P}}}/\sqrt{{b}^{2}+{\varepsilon }^{2}}$. The ${v}_{{\rm{P}}}^{-1}$ and K_0i factors in Equation (39) conspire to provide the near-resonance condition for maximum response,

$$\begin{eqnarray}&&n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi }\approx \displaystyle \frac{0.6\,{v}_{{\rm{P}}}}{\sqrt{{b}^{2}+{\varepsilon }^{2}}},\end{eqnarray} \tag{ 43 }$$

where b is the impact parameter of the encounter, given by Equation (42). From the top panels of Figure 7, it is clear that the peak response shifts to smaller v_P with increasing I_z . This is easy to understand from the fact that the corresponding vertical frequency, Ω_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, θ_P, but very mildly on the azimuthal angle, ϕ_P. Moreover, the middle panels indicate that more planar encounters (larger θ_P) induce stronger responses. This is a consequence of the fact that satellites spend more time near the disk during more planar encounters.

The in-plane structure of the disk response depends on the relative contribution of the different (l, m) modes. From Figure 8, 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). Because, for a Sgr-like encounter, f_1,nlm /f₀ ≳ 1 for these dominant modes, the response to the impact by Sgr is in fact nonlinear in the solar neighborhood. Either way, a satellite encounter is typically found to excite strong m = −2 modes, i.e., two-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.

Figure 9 elucidates that the bending mode response dominates for slower encounters, i.e., smaller v_P, and at guiding radii far from the impact radius, R_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_P encounters versus that of breathing modes for high-v_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 draining energy from those near the other turning point, thereby driving bending wave perturbations that are asymmetric about the midplane. 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 that 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_c ∼ 6–7 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 they 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 one or more dark satellites 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.

To evaluate the disk response to satellite encounters using Equation (21), we first evaluate the τ integral (with t_i → − ∞ ) of the satellite potential given in Equation (38) and then compute the Fourier transform of the result. This yields the expression for the response in Equation (21) with

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & = & \exp \left[-i{\rm{\Omega }}t\right]{\int }_{-\infty }^{t}d\tau \,\exp \left[i{\rm{\Omega }}\tau \right]\,{{\rm{\Phi }}}_{{nlm}}\left({\boldsymbol{I}},\tau \right)\\ & = & \displaystyle \frac{\exp \left[-i{\rm{\Omega }}t\right]}{{\left(2\pi \right)}^{3}}{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[-{{inw}}_{z}\right]{\int }_{0}^{2\pi }{{dw}}_{R}\exp \left[-{{ilw}}_{R}\right]{\int }_{0}^{2\pi }{{dw}}_{\phi }\exp \left[-{{imw}}_{\phi }\right]{\int }_{-\infty }^{t}d\tau \,\exp \left[i{\rm{\Omega }}\tau \right]\,{{\rm{\Phi }}}_{{\rm{P}}}(z,R,\phi ,\tau ),\end{array}\end{eqnarray} \tag{ D1 }$$

where

$$\begin{eqnarray}&&{\rm{\Omega }}=n{{\rm{\Omega }}}_{z}+l{{\rm{\Omega }}}_{R}+m{{\rm{\Omega }}}_{\phi }.\end{eqnarray} \tag{ D2 }$$

We perform the inner τ integral of Φ_P to obtain

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{-\infty }^{t}d\tau \,\exp \left[i{\rm{\Omega }}\tau \right]\,{{\rm{\Phi }}}_{{\rm{P}}}(z,R,\phi ,\tau ) & = & -\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[i\displaystyle \frac{{\rm{\Omega }}{ \mathcal S }}{{v}_{{\rm{P}}}}\right]{\int }_{-\infty }^{t-{ \mathcal S }/{v}_{{\rm{P}}}}d\tau \displaystyle \frac{\exp \left[i{\rm{\Omega }}\tau \right]}{\sqrt{{\tau }^{2}+\left({{ \mathcal R }}^{2}+{\varepsilon }^{2}\right)/{v}_{{\rm{P}}}^{2}}}\\ & = & -\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[i\displaystyle \frac{{\rm{\Omega }}{ \mathcal S }}{{v}_{{\rm{P}}}}\right]{\int }_{-\infty }^{\left({v}_{{\rm{P}}}t-{ \mathcal S }\right)/\sqrt{{{ \mathcal R }}^{2}+{\varepsilon }^{2}}}{dx}\displaystyle \frac{\exp \left[i\left({\rm{\Omega }}\sqrt{{{ \mathcal R }}^{2}+{\varepsilon }^{2}}/{v}_{{\rm{P}}}\right)x\right]}{\sqrt{{x}^{2}+1}}\\ & = & -\displaystyle \frac{2\,{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[i\displaystyle \frac{{\rm{\Omega }}{ \mathcal S }}{{v}_{{\rm{P}}}}\right]{K}_{0i}\left(\displaystyle \frac{{\rm{\Omega }}\sqrt{{{ \mathcal R }}^{2}+{\varepsilon }^{2}}}{{v}_{{\rm{P}}}},\displaystyle \frac{{v}_{{\rm{P}}}t-{ \mathcal S }}{\sqrt{{{ \mathcal R }}^{2}+{\varepsilon }^{2}}}\right).\end{array}\end{eqnarray} \tag{ D3 }$$

Here, K_0i is defined as

$$\begin{eqnarray}&&{K}_{0i}(\alpha ,\beta )=\displaystyle \frac{1}{2}{\int }_{-\infty }^{\beta }{dx}\,\displaystyle \frac{\exp \left[i\alpha x\right]}{\sqrt{{x}^{2}+1}},\end{eqnarray} \tag{ D4 }$$

which asymptotes to the zeroth order modified Bessel function of the second kind, ${K}_{0}\left(\left|\alpha \right|\right)$, in the limit β → ∞ . ${ \mathcal R }$ and ${ \mathcal S }$ are, respectively, the projections perpendicular and parallel to the direction of v _P of the vector connecting the point of observation, (z, R, ϕ), with the point of impact, and are given by

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}^{2} & = & {\left[R\sin (\phi -{\phi }_{{\rm{P}}})+{R}_{{\rm{d}}}\sin {\phi }_{{\rm{P}}}\right]}^{2}+{\left[(R\cos (\phi -{\phi }_{{\rm{P}}})-{R}_{{\rm{d}}}\cos {\phi }_{{\rm{P}}})\cos {\theta }_{{\rm{P}}}-z\sin {\theta }_{{\rm{P}}}\right]}^{2}\\ { \mathcal S } & = & (R\cos (\phi -{\phi }_{{\rm{P}}})-{R}_{{\rm{d}}}\cos {\phi }_{{\rm{P}}})\sin {\theta }_{{\rm{P}}}+z\cos {\theta }_{{\rm{P}}}.\end{array}\end{eqnarray} \tag{ D5 }$$

In deriving Equation (D3), we have only considered the direct term in the expression for Φ_P given in Equation (38); the indirect term turns out to be subdominant.

In the large time limit, i.e., $t\gg { \mathcal S }/{v}_{{\rm{P}}}$, K_0i asymptotes to ${K}_{0}\left(\left|{\rm{\Omega }}\right|\sqrt{{{ \mathcal R }}^{2}+{\varepsilon }^{2}}/{v}_{{\rm{P}}}\right)$. We substitute the expressions for R and z in terms of ( w , I ) given in Equations (C2) and (C4) in the above expressions for ${ \mathcal R }$ and ${ \mathcal S }$. Further substituting the resultant τ integral from Equation (D3) in Equation (D1), substituting w_ϕ in terms of ϕ using Equation (C3), adopting the small I_R limit, and performing the w_R integral, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & \approx & -\displaystyle \frac{2{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[-i{\rm{\Omega }}t\right]\times \exp \left[-i\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}\cos {\phi }_{{\rm{P}}}}{{v}_{{\rm{P}}}}{R}_{{\rm{d}}}\right]\times \exp \left[{il}{\tan }^{-1}\displaystyle \frac{2m{{\rm{\Omega }}}_{\phi }}{{R}_{{\rm{c}}}\kappa }\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\right]\\ & & \times \displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[-{{inw}}_{z}\right]\exp \left[i\displaystyle \frac{{\rm{\Omega }}\cos {\theta }_{{\rm{P}}}}{{v}_{{\rm{P}}}}z\right]{\int }_{0}^{2\pi }d\phi \,\exp \left[-{im}\phi \right]\,\exp \left[i\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}\cos \left(\phi -{\phi }_{{\rm{P}}}\right)}{{v}_{{\rm{P}}}}{R}_{{\rm{c}}}\right]\\ & & \times {J}_{l}\left(\sqrt{{\left(\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}}{{v}_{{\rm{P}}}}\right)}^{2}{\cos }^{2}\left(\phi -{\phi }_{{\rm{P}}}\right)+{\left(\displaystyle \frac{2m{{\rm{\Omega }}}_{\phi }}{{R}_{{\rm{c}}}\kappa }\right)}^{2}}\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\right){K}_{0i}\left(\displaystyle \frac{{\rm{\Omega }}\sqrt{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}}{{v}_{{\rm{P}}}},\displaystyle \frac{{v}_{{\rm{P}}}t-{{ \mathcal S }}_{{\rm{c}}}}{\sqrt{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}}\right),\end{array}\end{eqnarray} \tag{ D6 }$$

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

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

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t)\mathop{\approx }\limits^{{\rm{t}}\to \infty }-\displaystyle \frac{2{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[-i{\rm{\Omega }}t\right]\times \exp \left[-i\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}\cos {\phi }_{{\rm{P}}}}{{v}_{{\rm{P}}}}{R}_{{\rm{d}}}\right]\\ \,\times \,\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[-{{inw}}_{z}\right]\exp \left[i\displaystyle \frac{{\rm{\Omega }}\cos {\theta }_{{\rm{P}}}}{{v}_{{\rm{P}}}}z\right]{\int }_{0}^{2\pi }d\phi \,\exp \left[-{im}\phi \right]\,\exp \left[i\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}\cos \left(\phi -{\phi }_{{\rm{P}}}\right)}{{v}_{{\rm{P}}}}{R}_{{\rm{c}}}\right]\\ \,\times \,\exp \left[{il}{\tan }^{-1}\displaystyle \frac{2m{{\rm{\Omega }}}_{\phi }}{{R}_{{\rm{c}}}\kappa }\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\right]\\ \,\times \,\left[{\zeta }^{\left(0\right)}{J}_{l}\left(\chi \right)-i{\zeta }^{\left(1\right)}{J}_{l}^{{\rm{{\prime} }}}\left(\chi \right)-\displaystyle \frac{1}{2}{\zeta }^{\left(2\right)}{J}_{l}^{{\rm{{\prime} }}{\rm{{\prime} }}}\left(\chi \right)\right],\end{array}\end{eqnarray} \tag{ D7 }$$

where

$$\begin{eqnarray}&&\chi =\sqrt{{\left(\displaystyle \frac{{\rm{\Omega }}\sin {\theta }_{{\rm{P}}}}{{v}_{{\rm{P}}}}\right)}^{2}{\cos }^{2}\left(\phi -{\phi }_{{\rm{P}}}\right)+{\left(\displaystyle \frac{2m{{\rm{\Omega }}}_{\phi }}{{R}_{{\rm{c}}}\kappa }\right)}^{2}}\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }},\end{eqnarray} \tag{ D8 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{\zeta }^{(0)} & = & {K}_{0}\left(\eta \right),\\ {\zeta }^{(1)} & = & \sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\,\displaystyle \frac{\partial {{ \mathcal R }}_{{\rm{c}}}}{\partial {R}_{{\rm{c}}}}\displaystyle \frac{{{ \mathcal R }}_{{\rm{c}}}}{\sqrt{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}}\displaystyle \frac{\left|{\rm{\Omega }}\right|}{{v}_{{\rm{P}}}}{K}_{0}^{{\prime} }\left(\eta \right),\\ {\zeta }^{(2)} & = & \displaystyle \frac{2{I}_{R}}{\kappa }\left[{\left(\displaystyle \frac{\partial {{ \mathcal R }}_{{\rm{c}}}}{\partial {R}_{{\rm{c}}}}\right)}^{2}\displaystyle \frac{{{ \mathcal R }}_{{\rm{c}}}^{2}}{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}\displaystyle \frac{{{\rm{\Omega }}}^{2}}{{v}_{{\rm{P}}}^{2}}{K}_{0}^{{\prime\prime} }(\eta )+\left\{\displaystyle \frac{{\partial }^{2}{{ \mathcal R }}_{{\rm{c}}}}{\partial {R}_{c}^{2}}\displaystyle \frac{{{ \mathcal R }}_{{\rm{c}}}}{\sqrt{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}}+{\left(\displaystyle \frac{\partial {{ \mathcal R }}_{{\rm{c}}}}{\partial {R}_{{\rm{c}}}}\right)}^{2}\displaystyle \frac{{\varepsilon }^{2}}{{\left({{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}\right)}^{3/2}}\right\}\displaystyle \frac{\left|{\rm{\Omega }}\right|}{{v}_{{\rm{P}}}}{K}_{0}^{{\prime} }(\eta )\right],\end{array}\end{eqnarray} \tag{ D9 }$$

with

$$\begin{eqnarray}&&\eta =\displaystyle \frac{\left|{\rm{\Omega }}\right|\sqrt{{{ \mathcal R }}_{{\rm{c}}}^{2}+{\varepsilon }^{2}}}{{v}_{{\rm{P}}}}.\end{eqnarray} \tag{ D10 }$$

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

Substituting the expression for ${{ \mathcal I }}_{{nlm}}$ given in Equation (D7) in the expression for the disk response given in Equation (21), 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,⊙ = h_z σ_z,⊙ and R_c (L_z ) = 8 kpc and marginalized over I_R , are summarized in Table 1.

The disk response in the general case, expressed by Equation (39), depends on several encounter parameters (R_d, θ_P, and ϕ_P), and it 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_P with a Plummer sphere of mass M_P and size ε is given by

$$\begin{eqnarray}&&{\rm{\Delta }}E=4\pi {\left(\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{v}_{{\rm{P}}}}\right)}^{2}{\int }_{0}^{\infty }{I}_{0}^{2}(R){\rm{\Sigma }}(R)\displaystyle \frac{{dR}}{R},\end{eqnarray} \tag{ D11 }$$

where

$$\begin{eqnarray}&&{I}_{0}(R)={\int }_{1}^{\infty }\displaystyle \frac{{M}_{{\rm{P}}}(\zeta R)}{{M}_{{\rm{p}}}}\,\displaystyle \frac{d\zeta }{{\zeta }^{2}{\left({\zeta }^{2}-1\right)}^{1/2}}.\end{eqnarray} \tag{ D12 }$$

Knowing that the enclosed mass profile of a Plummer sphere is given by ${M}_{{\rm{P}}}(R)={M}_{{\rm{P}}}{R}^{3}{\left({R}^{2}+{\varepsilon }^{2}\right)}^{-3/2}$, we find that I₀(R) = R²/(R² + ε²), which yields

$$\begin{eqnarray}&&{\rm{\Delta }}E=4\pi {\left(\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{v}_{{\rm{P}}}}\right)}^{2}{\int }_{0}^{\infty }{\rm{\Sigma }}(R)\displaystyle \frac{{R}^{3}{dR}}{{\left({R}^{2}+{\varepsilon }^{2}\right)}^{2}}.\end{eqnarray} \tag{ D13 }$$

Now we compute the disk response to the face-on satellite encounter using Equations (21) and (D7)–(D10). For a perpendicular face-on impact through the center of the disk, we have R_d = 0 and θ_P = 0, implying that ${{ \mathcal R }}_{{\rm{c}}}$ becomes R_c. The corresponding response is greatly simplified. In the large time and small I_R limit, it is given by Equation (21) with

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & \approx & -\displaystyle \frac{2{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[-i{\rm{\Omega }}t\right]\,{\delta }_{m,0}\times \displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[-{{inw}}_{z}\right]\exp \left[i\displaystyle \frac{{\rm{\Omega }}z}{{v}_{{\rm{P}}}}\right]\\ & & \times \displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{{dw}}_{R}\,\exp \left[-{{ilw}}_{R}\right]{K}_{0}\left[\displaystyle \frac{\left|{\rm{\Omega }}\right|}{{v}_{{\rm{P}}}}\sqrt{{\varepsilon }^{2}+{\left({R}_{{\rm{c}}}+\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\sin {w}_{R}\right)}^{2}}\right],\end{array}\end{eqnarray} \tag{ D14 }$$

where the ϕ 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 ${{ \mathcal I }}_{{nlm}}$:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & \approx & i\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\exp \left[-i{\rm{\Omega }}t\right]\,{\delta }_{m,0}\left({\delta }_{l,1}-{\delta }_{l,-1}\right)\times \displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[-{{inw}}_{z}\right]\exp \left[i\displaystyle \frac{{\rm{\Omega }}z}{{v}_{{\rm{P}}}}\right]\\ & & \times \sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\displaystyle \frac{{R}_{{\rm{c}}}}{\sqrt{{\varepsilon }^{2}+{R}_{c}^{2}}}\,\displaystyle \frac{\left|{\rm{\Omega }}\right|}{{v}_{{\rm{P}}}}{K}_{0}^{{\prime} }\left[\displaystyle \frac{\left|{\rm{\Omega }}\right|}{{v}_{{\rm{P}}}}\sqrt{{\varepsilon }^{2}+{R}_{c}^{2}}\right].\end{array}\end{eqnarray} \tag{ D15 }$$

In the impulsive limit, v_P → ∞ , this becomes

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal I }}_{{nlm}}({\boldsymbol{I}},t) & \approx & i\,{\delta }_{n,0}{\delta }_{m,0}\left({\delta }_{l,1}-{\delta }_{l,-1}\right)\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\\ & & \times \,\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\displaystyle \frac{{R}_{{\rm{c}}}}{{\varepsilon }^{2}+{R}_{c}^{2}}\exp \left[-i\,l\kappa \,t\right],\end{array}\end{eqnarray} \tag{ D16 }$$

which can be substituted in Equation (21) to yield

$$\begin{eqnarray}&&{f}_{1{nlm}}\left({\boldsymbol{I}},t\right)={f}_{0}({\boldsymbol{I}})\times {\delta }_{n,0}{\delta }_{m,0}\left({\delta }_{l,1}-{\delta }_{l,-1}\right)\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\displaystyle \frac{l\kappa }{{\sigma }_{R}^{2}}\sqrt{\displaystyle \frac{2{I}_{R}}{\kappa }}\displaystyle \frac{{R}_{{\rm{c}}}}{{\varepsilon }^{2}+{R}_{c}^{2}}\exp \left[-i\,l\kappa \,t\right],\end{eqnarray} \tag{ D17 }$$

with f₀ given by Equation (15). Hence, the response is given by

$$\begin{eqnarray}\begin{array}{rcl}{f}_{1}\left({\boldsymbol{w}},{\boldsymbol{I}},t\right) & = & \displaystyle \sum _{n=-\infty }^{\infty }\displaystyle \sum _{l=-\infty }^{\infty }\displaystyle \sum _{m=-\infty }^{\infty }\exp \left[i({{nw}}_{z}+{{lw}}_{R}+{{mw}}_{\phi })\right]{f}_{1{nlm}}({\boldsymbol{I}},t)\\ & = & {f}_{0}({\boldsymbol{I}})\times \displaystyle \frac{2{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\displaystyle \frac{\sqrt{2\kappa {I}_{R}}}{{\sigma }_{R}^{2}}\displaystyle \frac{{R}_{{\rm{c}}}}{{\varepsilon }^{2}+{R}_{c}^{2}}\cos \left({w}_{R}-\kappa t\right),\end{array}\end{eqnarray} \tag{ D18 }$$

which shows that the satellite passage introduces a relative overdensity, ${f}_{1}\left({\boldsymbol{w}},{\boldsymbol{I}},t\right)/{f}_{0}({\boldsymbol{I}})$, that scales as $\sim {R}_{{\rm{c}}}/\left({\varepsilon }^{2}+{R}_{c}^{2}\right)$, which increases from zero at the center, peaks at R_c = ε, and asymptotes to zero again at large R_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 dE/dt = ∂E/∂ I · d I /dt, where ∂E/∂ I = Ω = (Ω_z , Ω_R , Ω_ϕ ) and d I /dt = − ∂Φ_P/∂ w from Hamilton's equations of motion. Thus, the total phase-averaged energy injected per unit of phase space can be obtained as follows:

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\left({\boldsymbol{I}}\right)\right\rangle =\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}\int d{\boldsymbol{w}}{\int }_{-\infty }^{\infty }{dt}\,\displaystyle \frac{{dE}}{{dt}}{f}_{1}({\boldsymbol{I}},t)=-\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}\int d{\boldsymbol{w}}{\int }_{-\infty }^{\infty }{dt}\,{\boldsymbol{\Omega }}\cdot \displaystyle \frac{\partial {{\rm{\Phi }}}_{{\rm{P}}}}{\partial {\boldsymbol{w}}}{f}_{1}({\boldsymbol{I}},t).\end{eqnarray} \tag{ D19 }$$

We can substitute the Fourier series expansions of Φ_P and f₁ given in Equations (10) in the above expression and integrate over w to obtain (Weinberg 1994a, 1994b)

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\left({\boldsymbol{I}}\right)\right\rangle =i\displaystyle \sum _{{nlm}}\left(n{{\rm{\Omega }}}_{z}+l\kappa +m{{\rm{\Omega }}}_{\phi }\right){\int }_{-\infty }^{\infty }{dt}\,{{\rm{\Phi }}}_{{nlm}}^{* }({\boldsymbol{I}},t){f}_{1{nlm}}({\boldsymbol{I}},t).\end{eqnarray} \tag{ D20 }$$

We can now substitute the form of Φ_P for a Plummer perturber given in Equation (38), with r _P and r given by Equations (36) and (37). The time integral can thus be written as

$$\begin{eqnarray}\begin{array}{rcl}{\int }_{-\infty }^{\infty }{dt}\,{{\rm{\Phi }}}_{{nlm}}^{* }({\boldsymbol{I}},t){f}_{1{nlm}}({\boldsymbol{I}},t) & = & -\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}{\int }_{0}^{2\pi }{{dw}}_{z}\exp \left[{{inw}}_{z}\right]{\int }_{0}^{2\pi }{{dw}}_{R}\exp \left[{{ilw}}_{R}\right]{\int }_{0}^{2\pi }{{dw}}_{\phi }\exp \left[{{imw}}_{\phi }\right]\\ & & \times {\int }_{-\infty }^{\infty }{dt}\displaystyle \frac{{{GM}}_{{\rm{P}}}}{\sqrt{{\left({v}_{{\rm{P}}}t-z\right)}^{2}+{R}^{2}+{\varepsilon }^{2}}}{f}_{1{nlm}}({\boldsymbol{I}},t).\end{array}\end{eqnarray} \tag{ D21 }$$

Using Equations (C2) and (C4) to express R and z in terms of ( w , I ), and substituting the form for f_1nlm ( I , t) from Equation (D17), we can perform the above integrals over w and t. Substituting the result in Equation (D20), we obtain

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\left({\boldsymbol{I}}\right)\right\rangle ={\left(\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\right)}^{2}{f}_{0}({\boldsymbol{I}})\displaystyle \frac{2\kappa {I}_{R}}{{\sigma }_{R}^{2}}\,\displaystyle \frac{{R}_{c}^{2}}{{\left({\varepsilon }^{2}+{R}_{c}^{2}\right)}^{2}}.\end{eqnarray} \tag{ D22 }$$

The total energy, ΔE_tot, imparted into the disk by the impulsive satellite passage can be computed by integrating the above expression over I and w (which simply introduces a factor of ${\left(2\pi \right)}^{3}$, because $\left\langle {\rm{\Delta }}E\left({\boldsymbol{I}}\right)\right\rangle$ is already phase-averaged), using Equation (15) and transforming from L_z to R_c using the Jacobian dL_z /dR_c = R_c κ²/2Ω_ϕ . This yields

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{tot}}=4\pi {\left(\displaystyle \frac{{{GM}}_{{\rm{P}}}}{{v}_{{\rm{P}}}}\right)}^{2}{\int }_{0}^{\infty }{{dR}}_{{\rm{c}}}{R}_{{\rm{c}}}{\rm{\Sigma }}({R}_{{\rm{c}}})\displaystyle \frac{{R}_{{\rm{c}}}^{2}}{{\left({\varepsilon }^{2}+{R}_{{\rm{c}}}^{2}\right)}^{2}}.\end{eqnarray} \tag{ D23 }$$

This is indeed the expression for ΔE_tot derived under the impulse approximation given by Equation (D13).