© The Authors 2025

1 Introduction

In the Gaia Data Release 2 (DR2; Gaia Collaboration 2016, 2018a) data, Antoja et al. (2018) discovered a ‘one-arm’ spiral pattern within the stellar distribution in the vertical position (z) versus vertical velocity (v_z) space. This discovery provides evidence of incomplete phase mixing, indicating that the Milky Way (MW) was perturbed relatively recently (≲ 1 Gyr) and has not yet reached equilibrium. The phase spiral has been studied extensively because it contains two kinds of important information about the Galactic dynamics: the nature of past perturbation event(s) (e.g. Binney & Schönrich 2018; Khoperskov et al. 2019; Laporte et al. 2019; Bland-Hawthorn et al. 2019) and the Galactic potential (e.g. Widmark et al. 2021a, b, 2022a, b; Guo et al. 2024).

While various mechanisms have been proposed to explain the origin of the one-arm phase spiral, the external perturbation by the Sagittarius dwarf (Sgr) is considered a plausible scenario and has been widely studied. Theoretical studies using analytical models (e.g. Binney & Schönrich 2018; Banik et al. 2022, 2023) and controlled N-body simulations (e.g. Laporte et al. 2018, 2019; Bland-Hawthorn & Tepper-García 2021; Hunt et al. 2022; Bennett et al. 2022) have confirmed the emergence of the phase spiral in galactic discs perturbed by Sgr-like satellites. Not only present in such controlled simulations, phase spirals have also been discovered in cosmological simulations (García-Conde et al. 2022; Grand et al. 2023). Grand et al. (2023) demonstrated that the dark-matter (DM) wake generated by a satellite, rather than its direct impact, can excite the phase spiral. Similarly, Laporte et al. (2018, 2019) showed that the torque from the DM wake is larger than that from the main body of the Sgr during the early stage of the MW-Sgr interaction. Theoretical studies further suggest that Sgr induces other phase-space substructures and the Galactic warp (e.g. Gómez et al. 2013; Khanna et al. 2019; Tepper-García et al. 2022; Binney 2024). Apart from the Sgr scenario, other phase spiral formation mechanisms have been proposed, including bar buckling (Khoperskov et al. 2019) and the combined effects of many small perturbations (Tremaine et al. 2023), for instance, due to DM sub-halos (Gilman et al. 2025). However, their contributions to other distorted structures such as the warp are unclear.

Hunt et al. (2022) discovered another type of phase spiral in the Gaia DR3 data (Gaia Collaboration 2023), finding that stars with small guiding radii exhibit the ‘two-arm’ phase spiral. The one-arm and two-arm phase spirals correspond to different global vertical oscillation modes in the Galactic disc. The one-arm phase spiral is associated with the ‘bending’ mode, where the regions above and below the mid-plane oscillate in the same direction. In contrast, the two-arm phase spiral is associated with the ‘breathing’ mode, where the regions above and below the mid-plane oscillate in opposite directions. Although the one-arm phase spiral is observed in a wide area of the Galactic disc (e.g. Wang et al. 2019; Xu et al. 2020, 2023; Antoja et al. 2023), the two-arm phase spiral is detected only in specific radial and angular momentum ranges (Alinder et al. 2024); therefore, two-arm phase spirals are generally considered to arise from local internal processes rather than global external perturbations (Candlish et al. 2014; Hunt et al. 2022; Li et al. 2023; Chiba et al. 2025). Among potential internal drivers, spiral arms are plausible origins of the two-arm phase spiral, as they naturally excite the breathing mode (Debattista 2014; Faure et al. 2014; Monari et al. 2016a, b; Khachaturyants et al. 2022; Ghosh et al. 2022; Kumar et al. 2022; Asano et al. 2024).

Beyond the phase spirals, Gaia and other surveys have provided further evidence of bending and breathing oscillations on both local and global scales (e.g. Widrow et al. 2012; Williams et al. 2013; Xu et al. 2015; Gaia Collaboration 2018b; Kawata et al. 2018; Schönrich & Dehnen 2018; Wang et al. 2018a, b, 2020; Carrillo et al. 2019; Bennett & Bovy 2019; Khanna et al. 2019; López-Corredoira et al. 2020; Cheng et al. 2020; McMillan et al. 2022; Ghosh et al. 2022; Widmark et al. 2022c; Laporte et al. 2022; Ardèvol et al. 2023; Asano et al. 2024; Wang et al. 2024; Lin et al. 2024; Poggio et al. 2025). Many studies have investigated internal and external processes separately, but in reality the vertical oscillations in the MW disc should be driven by a combination of both. Bennett & Bovy (2021) suggested that the vertical asymmetry in the solar neighbourhood cannot be fully explained by the Sgr impact alone. Moreover, internal and external processes are not entirely independent, as external perturbations influence the internal structures of the disc. For instance, tidal interactions with satellite galaxies induce spiral arms in host discs (Antoja et al. 2022), which can further perturb the motion of disc stars. To fully understand the excitation and evolution of the bending mode and the breathing mode, it is crucial to consider the interplay between internal and external processes.

N-body simulations are useful tools to study the self-consistent evolution of galactic discs. In the context of the Sgr impact, most of the previous studies used dynamically hot discs, where spontaneous bar and spiral formations are suppressed to isolate the effects of external perturbation (but see Bland-Hawthorn & Tepper-García 2021). However, such hot disc models underestimate the influence of self-gravity, which affects the amplitude and winding rate of the phase spiral (Darling & Widrow 2019a, b; Widrow 2023). To address these limitations and bridge the gap between theory and observations, we employ a ‘cold’ disc model that accounts for the influence of self-gravity and better replicates the dynamical conditions of the MW. Our high-resolution simulations, powered by a graphics processing unit (GPU)-based code (Bédorf et al. 2012, 2014), utilise five billion particles to resolve fine phase-space structures and directly compare the results with Gaia data. Unlike previous billion-particle N-body models (Hunt et al. 2021; Bennett et al. 2022; Stelea et al. 2024), our disc model is colder and incorporates disc and bulge properties more consistent with the real MW. This approach allows us to investigate the disc’s dynamical evolution under external perturbations caused by a Sgr-like satellite galaxy while also capturing the potential role of internal processes such as spiral arms.

This paper is organised as follows. In Section 2, we describe the details of the N-body simulations. In Section 3, we analyse the global vertical oscillations of the MW disc induced by repeated impacts from the satellite, focusing on the bending and breathing modes and their time evolution. In Section 4, we isolate the effects of a single impact to clarify the role of the dwarf galaxy’s perturbation in driving the global vertical oscillations, by examining the long-term evolution of the bending and breathing modes after a single satellite passage. In Section 5, we shift to the local scale, investigating the spatial variation and time evolution of phase spirals, which link global vertical oscillations to local phase-space structures. In Section 6, we discuss what we can infer about the Galactic disc structure and its evolution history from the real observation results and our simulation. Finally, Section 7 provides a summary of our findings.

2 N-body simulation

2.1 Initial condition and simulation

We performed N-body simulations of a MW-like disc galaxy perturbed by a Sgr-like satellite galaxy. The MW model is identical to MWa of Fujii et al. (2019). We generated the initial condition from the same parameters as MWa using Galactics (Kuijken & Dubinski 1995; Widrow & Dubinski 2005; Widrow et al. 2008). It comprises a DM halo, a classical bulge, and a disc.

The DM halo follows the Navarro-Frenk-White (NFW; Navarro et al. 1997) profile:

$${\rho }_{\mathrm{h}}(r)=\frac{{\rho }_{\mathrm{h}0}}{\left(r/a_{\mathrm{h}}\right){(1+r/a_{\mathrm{h}})}^2},$$(1)

where a_h and ρ_h0 are the scale radius and the characteristic density, respectively. Here, a_h is set to 10 kpc, and ρ_{h, 0} is determined through the characteristic velocity dispersion ${\sigma }_{\mathrm{h}}={\left(4\pi G{a}_{\mathrm{h}}^2{\rho }_{\mathrm{h}0}\right)}^{1/2}=420\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$, where G is the gravitational constant. The truncation parameter (Widrow & Dubinski 2005) is set to ɛ_h = 0.85. The rotation fraction, α_h, dictates the spin of the halo. This parameter represents the fraction of particles rotating in the same direction as the disc. Thus, if α_h equals 0.5, the halo exhibits no rotation. The chosen value of α_h = 0.8 results in a pattern speed and bar length consistent with the observed values (Bland-Hawthorn & Gerhard 2016).

The classical bulge follows the Hernquist profile (Hernquist 1990):

$${\rho }_{\mathrm{b}}(r)=\frac{{\rho }_{\mathrm{b}0}}{\left(r/a_{\mathrm{b}}\right){(1+r/a_{\mathrm{b}})}^3}.$$(2)

The scale radius, a_b, the characteristic velocity dispersion, ${\sigma }_{\mathrm{b}}={\left(4\pi G{a}_{\mathrm{b}}^2{\rho }_{\mathrm{b}0}\right)}^{1/2}$, and the truncation parameter, ɛ_b, are set to 0.75 kpc, 330 km s⁻¹, and 0.99, respectively. The rotation fraction, α_b, is set to 0.5, indicating that the classical bulge does not initially rotate.

The density distribution of the disc is given by a radially exponential and vertically isothermal profile:

$${\rho }_{\mathrm{d}}(R,z)={\rho }_{\mathrm{d}0}\exp (-\frac{R}{R_{\mathrm{d}}}){\mathrm{sech}}^2\left(\frac{z}{z_{\mathrm{d}}}\right).$$(3)

The five parameters: scale radius, R_d, the scale height, z_d, the total mass, M_d, the disc truncation radius, R_out, and the sharpness of truncation, δR_out, characterise the model. They are set to R_d = 2.3 kpc, z_d = 0.2 kpc, M_d = 3.61 × 10¹⁰ M_⊙, R_out = 30 kpc, and δ R_out = 0.8 kpc, respectively. The radial velocity dispersion follows ${\sigma }_R(R{)}^2={\sigma }_{R0}^2\exp (-R/R_{\mathrm{d}})$, with σ_R0 = 94 km s⁻¹.

The particle numbers for the DM halo, bulge, and disc components are N_h = 4.9 B, N_b = 3.1 M, and N_d = 213 M, respectively. We evolved the MW model for 8 Gyr and used the final snapshot for the initial condition of the host galaxy in the perturbed simulation. The initial condition parameters for the host galaxy are summarised in Table 1.

The dwarf model comprises both DM and stellar components. We generated the initial condition with Agama (Vasiliev 2019). The distribution of the DM particles follows the NFW profile. The total mass and the scale radius of the DM halo are 5 × 10¹⁰ M_⊙ and 7.5 kpc, respectively. The outer cut-off radius is set to 33 kpc, which corresponds to the tidal radius of the dwarf at the initial position. The stellar component follows the Hernquist profile. The total mass and the scale radius are 1 × 10⁹ M_⊙ and 2 kpc, respectively. Neither the DM component nor the stellar component has a spin. The numbers of the DM and stellar particles are N_DM = 297 M and N_star = 5.8 M, respectively. The initial condition parameters for the dwarf galaxy are summarised in Table 2.

Combining the MW-like host with the Sgr-like satellite, we simulated their interaction for ∼ 3 Gyr. We used the parallel GPU tree-code Bonsai¹ (Bédorf et al. 2012, 2014) on Pegasus at Center for Computational Sciences, University of Tsukuba. We adopted a softening length of 0.01 kpc, an opening angle of 0.4, and a shared integration timestep of 0.61 Myr. Snapshots were output at intervals of 9.78 Myr. The simulation took about 10 hours with 43 NVIDIA H100 GPUs. The simulation data is available at http://galaxies.astron.s.u-tokyo.ac.jp.

2.2 The time evolution of the dwarf

Fig. 1 shows the orbit and mass evolution of the dwarf. The dwarf is initially located at the apocentre of r ≈ 130 kpc and completes three orbit loops by the end of the simulation, with its orbit decaying due to dynamical friction. The dwarf’s position at t = 1.78 Gyr is closest to the present-day position of the Sgr (Vasiliev & Belokurov 2020). The vertical red line in the figure indicates this time. The lower panel shows the time evolution of the dwarf’s mass. The solid and dashed lines indicate the total mass (M_total=M_DM+M_star) and the stellar mass (M_star) within 5 kpc from the centre of the dwarf, respectively. The centre is defined as the density peak of the stellar component. The dotted line represents the stellar mass fraction with respect to the total mass (M_star/M_total). It exhibits staircase patterns, reflecting the efficient tidal stripping at the pericentre. In the final epoch of the simulation, the total mass is 8 × 10⁸ M_⊙, which is comparable with the observationally estimated Sgr mass, ∼ 5 × 10⁸ M_⊙ (Vasiliev & Belokurov 2020). Fig. 2 presents the projections of the dwarf’s orbital trajectory into the x–y and x–z planes. The dwarf follows a nearly polar orbit similar to that of the Sgr. The cyan point and the cyan arrow indicate the position and velocity of the dwarf at t = 1.78 Gyr, respectively. The red square and the red dashed arrow indicate those of the Sgr (Vasiliev & Belokurov 2020). The purple dots represent stellar particles stripped from the dwarf, forming a stellar stream which roughly follows the orbit of the dwarf. This feature is qualitatively consistent with the observed Sgr stream (e.g. Majewski et al. 2003) as well as results of other simulations (e.g. Wang et al. 2022).

Fig. 1 Upper panel: position of the dwarf as a function of time. Blue solid, orange dashed, and green dotted lines show $r=\sqrt{x^2+y^2+z^2},R=\sqrt{x^2+y^2}$, and z, respectively. Lower panel: mass of the dwarf as a function of time. Blue solid and orange dashed lines show the total mass and stellar mass enclosed within 5 kpc. The Green dotted line shows the stellar mass fraction. The red vertical line indicates the time of t = 1.78 Gyr when the dwarf’s position is close to the present-day position of the Sgr. The final total mass is consistent with the presentday Sgr mass.

3 Global vertical oscillations from repeated impacts

In this section, we investigate the global vertical oscillations of the MW disc induced by repeated impacts from the Sgr-like dwarf galaxy. We first provide a qualitative description of the face-on evolution of the disc structure and then perform a quantitative analysis of the bending and breathing modes, focusing on their time evolution.

3.1 Face-on maps of the galactic disc

The external perturbation by the dwarf causes a shift in the centre of mass and a tilt of the disc with respect to the initial galactic plane. To focus on the internal structure of the disc, we redefined the coordinate system by setting (x, y, z) = (0, 0, 0) at the centre of mass of the host’s stellar components (bulge and disc) and aligning the z-axis with their total angular momentum vector.

Fig. 3 shows face-on maps of the galactic disc at t= 1.15 Gyr, which is ∼ 250 Myr after the first pericentre passage of the dwarf². The upper left panel is colour-coded by the particle numbers in 0.2 × 0.2 kpc² bins, which represent the surface density of the disc. The other panels show the mean azimuthal velocity (v_ϕ) with respect to the azimuthally averaged velocity $\left({\overline{v}}_{\varphi }\right)$, the mean radial velocity (v_R), the mean vertical coordinate (z), the mean vertical velocity (v_z), and the breathing velocity which is defined as the difference in the mean vertical velocities above and below the mid-plane (v_{z, > 0}–v_{z, < 0}). To reduce noise, the breathing velocity map is smoothed with a Gaussian filter of 0.2 kpc.

The surface density map (upper left panel) exhibits a bar and spiral arms. While the bar exists before the dwarf’s interaction, the prominent two-arm spiral pattern is induced by the perturbation from the dwarf. Additionally, we can identify faint arms. They are transient dynamic arms, which result from the disc’s self-instability not from the external perturbation.

The azimuthal and radial velocity maps (upper middle and upper right panels) exhibit quadrupole patterns in the central region of the galaxy, which are signatures of the bar. Beyond the bar radius, the two-arm patterns align with the tidally induced spiral arms. As studied by Antoja et al. (2022), a satellite fly-by initially induces a global quadrupole velocity pattern in the disc, which corresponds to the initial conditions for forming a two-arm kinematic density wave (Kalnajs 1973). This wave develops into the two-arm spiral pattern due to the differential rotation of the disc.

The vertical coordinate and the vertical velocity maps (lower left and lower middle panels) exhibit a corrugation of the galactic disc. The amplitude of the corrugation increases with the galactocentric radius. During the early stage of interaction, the galactic disc exhibits an S-shaped warp. It starts to emerge in the outer galaxy ∼ 200 Myr before the dwarf’s pericentre passage. The amplitude of the warp grows until the pericentre passage, after which it winds up rapidly, evolving into a corrugated structure whose amplitude decreases over time. This damping is faster in the inner galaxy (R ≲ 10 kpc) than in the outer galaxy (R ≳ 10 kpc).

In the breathing map (lower right panel), a two-arm pattern is observed along the tidally induced spiral arms. As demonstrated by previous theoretical studies (e.g. Debattista 2014), spiral arms naturally excite the breathing mode. Kumar et al. (2021, 2022) also performed N-body simulations of disc galaxies interacting with fly-by satellites and identified a prominent largescale breathing pattern excited by tidally influenced spiral arms. Figure 1 of Haines et al. (2019) shows similar face-on breathing maps from Laporte et al. (2018)’s N-body model, where we can identify two-arm breathing pattern possibly excited by tidally induced arms. In the central region, a quadrupole breathing pattern arises from the bar, while in the outermost region (R ≳ 15 kpc), a more complex breathing pattern emerges. This complexity likely results from the combination of the directly excited breathing mode due to the satellite impact and the spiralinduced breathing mode, as the former cannot be neglected in the outer galaxy (Banik et al. 2023).

Fig. 2 Projection of the dwarf’s orbit in the x–y plane (upper panel) and in the x–z plane (lower panel). The cyan line indicates the orbital trajectory of the dwarf. The cyan point and cyan solid arrow indicate the dwarf’s position and velocity at t = 1.78 Gyr, respectively. The red square and red dashed arrow indicate those for the Sgr at the present day estimated by Vasiliev & Belokurov (2020). Purple dots show the distribution of the stellar particles stripped away from the dwarf. The surface density of the host disc and the solar position are displayed on the background. The dwarf follows a nearly polar orbit similar to the Sgr, leaving a stellar stream of stripped particles along its orbital path.

Fig. 3 Face-on maps of the disc at t = 1.15 Gyr, which is ∼ 250 Myr after the first pericentre passage of the dwarf. Tidal forces from the dwarf induce prominent two-arm spiral arms and cause vertical corrugations (bending) and a breathing pattern in the galactic disc. The breathing pattern is aligned with the spiral arms.

3.2 Bending and breathing modes

We quantitatively evaluate the amplitudes of the bending and breathing modes and investigate their time evolution. Based on the symmetry of these two modes, we can evaluate their amplitude by comparing the vertical velocities above and below the mid-plane (Widrow et al. 2014; Khachaturyants et al. 2022; García-Conde et al. 2024). We defined the bending and breathing velocities as

$$V_{\text{bend}}(R,\varphi ,z,t)=\frac{1}{2}[{\overline{v}}_z(R,\varphi ,z,t)+{\overline{v}}_z(R,\varphi ,-z,t)],$$(4)

$$V_{\text{breath}}(R,\varphi ,z,t)=\frac{1}{2}[{\overline{v}}_z(R,\varphi ,z,t)-{\overline{v}}_z(R,\varphi ,-z,t)],$$(5)

where ${\overline{v}}_z(R,\varphi ,z,t)$ is the mean vertical velocity at position (R, ϕ, z) at time t. To obtain a vertically averaged value for each mode, V_n(n ∈ bend, breath), we calculated

$$⟨V_n⟩(R,\varphi ,t)=\frac{1}{Z_{\text{max}}}{\int }_0^{Z_{\text{max}}}{V}_n(R,\varphi ,z,t)\mathrm{d}z,$$(6)

where Z_max = 1.5 kpc. We then computed the azimuthally averaged squared amplitudes of the modes:

$$A_n^2(R,t)=\frac{1}{2\pi }{\int }_0^{2\pi }{⟨V_n⟩}^2(R,\varphi ,t)\mathrm{d}\varphi .$$(7)

For numerical evaluation, we applied binning with radial, azimuthal, and vertical resolutions of 0.5 kpc, 1°, and 0.1 kpc, respectively.

The first panel of Figure 4 shows the time evolution of the bending and breathing amplitudes (A_n) in the radial range R = 8–8.5 kpc. The vertical lines indicate the dwarf’s pericentre passages. Both the bending and breathing amplitudes begin to increase around the time of the first pericentre passage (t ∼ 0.9 Gyr), indicating that external perturbation triggers both modes. However, their responses differ in detail. The bending mode (solid line) increases sharply at t ∼ 0.85 Gyr (∼ 50 Myr before the pericentre passage) and decays shortly after the pericentre passage. In contrast, the breathing mode (dashed line) starts to increase after the pericentre passage, with its amplitude reaching saturation around t = 1.2 Gyr. Initially, the bending mode dominates, but approximately 300 Myr after the pericentre passage, the breathing mode becomes dominant. This transition of the dominant mode reflects the differences in the growth and decay timescales of the two modes, as discussed in more detail in the next section.

After the second pericentre passage (t ∼ 1.75 Gyr), the bending mode again shows a sharp increase in amplitude, whereas the breathing mode shows no clear response. During this epoch, the breathing amplitude gradually decreases. Between the second and third pericentres, 1.75 ≲ t ≲ 2.4 Gyr, the bending amplitude remains higher or comparable to the breathing amplitude. The third pericentre passage does not significantly impact either the bending or breathing modes.

Fig. 4 Time evolution of the bending mode, breathing mode, and surface density in the ring of R = 8–8.5 kpc. First panel: total amplitudes of the bending (solid line) and breathing (dashed line) as functions of time. Second panel: Fourier amplitudes of the bending mode for m = 0, 1, 2, 3, and 4. Different line styles correspond to different Fourier components as indicated in the legend. Third panel: same as the second panel but for the breathing mode. Fourth panel: time evolution of the Fourier amplitudes of the surface density. In all panels, vertical lines indicate the timing of the pericentre passages of the dwarf. The bending mode responds sharply to pericentre passages and decays quickly, while the breathing mode increases more slowly and persists longer, leading to a transition from bending-dominated to breathing-dominated state. The m = 1 and m = 2 are the dominant Fourier components for the bending and breathing modes, respectively.

3.3 Fourier decomposition

To gain deeper insight into the temporal evolution of the bending and breathing modes, we performed the Fourier decomposition of the bending and breathing modes. Similar Fourier analyses were employed in previous studies (e.g. Chequers & Widrow 2017; Chequers et al. 2018; Poggio et al. 2021; Grion Filho et al. 2021; García-Conde et al. 2024). The vertically averaged bending and breathing velocities of Equation (6) were decomposed as

$$⟨V_n⟩(R,\varphi ,t)=\sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}{Q}_{n,m}(R,t)\exp (im\varphi )$$(8)

$$=\sum _{m=0}^{\mathrm{\infty }}{A}_{n,m}(R,t)\cos \left\{m[\varphi -{\varphi }_{n,m}(R,t)]\right\}$$(9)

where Q_n,m, A_n,m, and ϕ_n,m are the complex coefficient, amplitude, and phase of the m-th Fourier component, respectively. Since ⟨V_n⟩ is a real function, the complex Fourier coefficients are related to the amplitude and the phase as

$$Q_{n,m}(R,t)=\{\begin{array}{ll}{A}_{n,0}(R,t)& \text{for} m=0,\\ \frac{1}{2}{A}_{n,m}(R,t)e^{-im{\varphi }_{n,m}(R,t)}& \text{for} m=1,2,\dots \\ \frac{1}{2}{A}_{n,-m}(R,t)e^{im{\varphi }_{n,-m}(R,t)}& \text{for} m=-1,-2,\dots \end{array}$$(10)

Similarly, we performed the Fourier decomposition of the normalised surface density,

$$\frac{\mathrm{\Sigma }(R,\varphi ,t)}{{\mathrm{\Sigma }}_0(R,t)}=\sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}{Q}_{\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s},m}(R,t)\exp (im\varphi )$$(11)

$$=1+\sum _{m=1}^{\mathrm{\infty }}{A}_{\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s},m}(R,t)\cos \left\{m[\varphi -{\varphi }_{\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s},m}(R,t)]\right\}$$(12)

where Σ₀(R, t) is the azimuthally averaged surface density.

The second panel of Fig. 4 shows the Fourier amplitudes of the bending mode for m = 0, 1, 2, 3, and 4. The m = 1 and m = 2 components correspond to the S-shaped and U-shaped warps, respectively. The m = 0 component represents the ring-shaped vertical oscillation. It is sometimes referred to as the bell mode (Merritt & Sellwood 1994; Sellwood & Merritt 1994). Among these, the m = 1 component is the most dominant, growing at t ∼ 0.85 Gyr, reaching its peak at t ∼ 0.95 Gyr, and subsequently undergoing a damping oscillation.

The m = 2 component, the second most dominant, also starts growing at t ∼ 0.85 Gyr. However, it lags behind the m = 1 component, reaching its maximum amplitude at t ∼ 1.2 Gyr. Interestingly, the m = 2 component experiences a secondary increase in amplitude after its initial decay just before the second pericentre passage of the dwarf, although the cause of this behaviour is unclear. In other similar simulations, the m = 2 component dominates over the m = 1 component. In the N-body simulation of Laporte et al. (2018), m = 2 component becomes dominant across a wide radial range 200−300 Myr after the first disc crossing (Poggio et al. 2021). Similarly, Bland-Hawthorn & Tepper-García (2021) show that the m = 2 component is the dominant bending component in the inner galaxy (R ≲ 10 kpc) in their simulation. The difference in the dominant Fourier component may be due to the differences in satellite orbits, since the relative strength of the m = 1 and m = 2 components depends on the satellite’s velocity (Banik et al. 2023).

After the second pericentre passage, the m = 1 component responds sensitively to external perturbation, while the amplitudes of the other components remain largely unchanged and comparable to each other. The third pericentre passage does not significantly affect any of the Fourier components.

The third panel of Fig. 4 illustrates the Fourier amplitudes of the breathing mode over time. The m = 2 component dominates, consistent with the total amplitude shown in the first panel. While the first pericentre passage triggers the m = 2 component, the second and third passages have little impact. The m = 2 breathing amplitude slowly decreases from t ∼ 1.7 Gyr and drops rapidly at t ∼ 2.3 Gyr.

The fourth panel of Fig. 4 shows the Fourier amplitudes of the surface density. The m = 2 component grows after the first pericentre passage, indicating the formation of tidally induced spiral arms, but its amplitude decreases around the second pericentre passage. It remains unclear whether this decline is directly caused by the second passage, as the trend begins earlier. A slight resurgence is observed after the third pericentre passage. While single fly-by interactions are known to generate two-arm spiral arms in smooth discs (Struck et al. 2011; Pettitt et al. 2016; Bland-Hawthorn & Tepper-García 2021; Antoja et al. 2022), it is uncertain whether such perturbations strengthen or weaken preexisting arms. Notably, the similarity in the Fourier amplitudes of the breathing mode and the surface density suggests that the breathing mode is excited by spiral arms.

3.4 Bending and breathing amplitudes as functions of time and radius

In the previous subsection, we examined the temporal evolution of the bending and breathing modes at a fixed radius. Here, we extend our analysis to explore how these modes vary as functions of both time and radius. This allows us to investigate the spatial dependence of their excitation and decay processes.

The first panel of Fig. 5 shows the bending amplitudes, as defined in Equation (7), as functions of time and the galactocentric radius. The first close encounter excites the bending mode across a wide radial range. The amplitude increases with R, as the dwarf crosses the outer disc (R = 25.5 kpc). The maps reveal ridge structures, indicating that the bending amplitude oscillates with time. These ridges gradually fade over time, suggesting that the bending oscillation behaves as a damping wave rather than a stationary wave. The second pericentre passage excites the bending mode again, but the ridge pattern becomes less regular compared to that after the first passage, reflecting increased complexity from repeated interactions. While the third pericentre passage causes a slight increase in the bending amplitude, its impact is notably weaker than the earlier passages.

The second panel of Fig. 5 shows the same plot but for the breathing mode. Before the dwarf’s first pericentre passage, a bar-induced breathing mode is visible as double peaks at R ∼ 1 kpc and R ∼ 2.5 kpc. These features are linked to the box y-peanut shape of the bar (Asano et al., in prep.). Outside the bar region, the breathing mode grows after the first pericentre passage, driven by tidally induced spiral arms, as observed in Figs. 3 and 4. Unlike the bending mode, the breathing mode exhibits slower growth and decay, persisting over longer timescales. For example, at R = 8 kpc, the breathing mode remains significant from t ∼ 1 Gyr to ∼ 2 Gyr. The second and third pericentre passages, however, do not strongly influence the breathing mode.

The third panel of Fig. 5 is colour-coded by the amplitude ratio between the bending and the breathing modes. In this map, red and blue regions represent dominance by the bending and breathing modes, respectively. Before the dwarf’s first pericentre passage, the breathing amplitude exceeds the bending amplitude across most parts of the disc, except in the outermost region. At the first pericentre passage, the bending mode becomes dominant. Over time, however, the dominant mode transitions back to the breathing mode, with this transition occurring faster in the inner galaxy than in the outer galaxy. By the time of the second pericentre passage, the breathing mode dominates within R ≲ 10 kpc.

This temporal transition reflects the distinct excitation and decay behaviours of the two modes. While the bending mode is excited instantaneously by the satellite’s perturbation and decays rapidly, the breathing mode is continuously driven by the spiral arms and persists as long as the arms remain sufficiently strong. This difference in temporal behaviour between the two modes drives the transition of the dominant mode.

Between the second and third pericentre passages, no clear transition is observed, since the breathing mode also begins to decay during this period. After the third pericentre passage, the outer galaxy (R ≳ 10 kpc) is dominated by the bending mode, while the inner galaxy (R ≲ 3 kpc) is dominated by the barinduced breathing mode. At intermediate radii, the amplitudes of the two modes are comparable, reflecting a balance between the decaying bending and breathing modes.

Fig. 5 Bending amplitude, breathing amplitude, and their ratio as functions of time and galactocentric radius. First panel: total amplitude of the bending mode. Second panel: total amplitude of the breathing mode. Third panel: amplitude ratio between the bending mode and the breathing mode. The vertical lines indicate the times of the dwarf’s pericentre passages. The dominant mode changes from bending to breathing. This transition occurs faster in the inner galaxy.

4 Long-term evolution of the bending mode and the breathing mode after a single impact

In the previous section, we investigated how repeated interactions with the dwarf galaxy introduce complexities into the MW disc. To isolate the impact of a single perturbation, we now investigate the long-term evolution of global vertical oscillations in the disc after a single impact. To this end, we performed the same N-body simulation of the MW-Sgr system but stopped the run at t = 1.2 Gyr, when the dwarf reaches its first apocentre. At this point, we removed the dwarf’s stellar and DM particles within 30 kpc from its centre. After removing the main body of the dwarf, we continued the simulation for an additional 1.5 Gyr.

4.1 Damping of the bending and breathing modes

The first panel of Fig. 6 shows the ratio of the bending to breathing amplitude as a function of t–t₀ and R in the single-impact model, where t₀ ∼ 0.9 Gyr is the pericentre passage time of the dwarf. Shortly after the pericentre passage, the bending mode dominates across the entire disc, but over time, the breathing mode becomes dominant, particularly in the inner regions (Fig. 5). The dashed curve represents R/σ_R(R), where σ_R(R) is the radial velocity dispersion at R. It roughly separates the bending-dominated and breathing-dominated regions. By analogy with wave damping in horizontally uniform slabs (Banik et al. 2022), we can interpret the timescale R/σ_R(R) as being relevant to horizontal mixing. We describe a mechanism in detail below.

The second panel of Fig. 6 shows the bending amplitude, normalised by its maximum value at each radius, as a function of t and R. The solid curves in the panel show (1, 3, 5, 7) × π/v(R)– 50 Myr, where ν(R) is the vertical frequency at R. To compute v, we constructed an axisymmetric potential model based on the parameters for the initial condition described in Section 2.1 using Agama (Vasiliev 2019). The bending amplitude shows a characteristic damping oscillation, with ridges of positive slope closely following the solid lines that represent 2π/ν(R). Over time, the amplitude decays by ∼ 40–60%, with the damping timescale corresponding to R/σ_R.

The damping observed here is consistent with horizontal (lateral)³ phase mixing in a horizontally uniform slab, where the response to an impulsive perturbation decays as

$$\exp (-\frac{{\sigma }^2{k}^2{t}^2}{2}),$$(13)

where σ and k are the horizontal velocity dispersion and the horizontal wave number, respectively, and the characteristic timescale of the damping is expressed as τ_D = 1 /(σ k) (Banik et al. 2022). Although our model is more complex than the infinite slab, the damping timescale due to horizontal mixing should be of the same order as that obtained by substituting σ ∼ σ_R and k ∼ 1/R into this formula (i.e. τ_D=R/σ_R).

Similar damping behaviour is observed in simpler systems. Section 5.2 of Binney & Tremaine (2008) examines the responses of an infinite, homogeneous stellar system with a Maxwellian distribution function. In this case, the spatial Fourier transform of the polarization function is proportional to Eq. (13). In general, waves in stellar systems are damped by two processes: phase mixing and Landau damping. Phase mixing is a kinematic process and occurs even in the absence of self-gravity, whereas Landau damping arises from self-gravitating effects and occurs more slowly than phase mixing. In the inner galaxy, where self-gravity is stronger, the damping timescale is longer than predicted by the R/σ_R scaling, as seen in the second panel of Fig. 6.

The third panel of Fig. 6 shows the breathing amplitude, normalised by its maximum values at each radius, as a function of t and R. The solid lines represent (1, 2, 3, 4) × π/ Ω (R), where Ω(R) is the circular frequency at R. The breathing mode takes approximately 1−1.5 circular periods to reach significant amplitudes after the pericentre passage. Once excited, it persists for up to ∼ 1 Gyr, driven by the tidally induced spiral arms responsible for exciting the mode. Unlike the bending mode, the breathing mode does not decay rapidly. As a result, the time evolution of the amplitude ratio between the two modes is primarily governed by the damping of the bending mode. While the bending mode fades, the breathing mode’s persistence reflects the longevity of the spiral arms. This distinction highlights the contrasting roles of phase mixing and spiral-arm dynamics in shaping the vertical oscillations of the disc.

Fig. 6 Evolution of the bending mode and the breathing mode in the single-impact model. First panel: bending and breathing amplitude ratio as a faction of time and the galactocentric radius. The horizontal axis indicates the time since the pericentre passage of the dwarf. A vertical dotted line indicates the pericentre time t–t0 = 0. Dashed curve represents R/σR(R). Second panel: bending amplitude as a function of time and the galactocentric radius. The amplitude is normalised by the maximum amplitude at each radius. Solid curves represent (1, 3, 5, 7) × π/v(R)−50 Myr. Third panel: same as the second panel but for the breathing mode. Solid lines represent (1, 2, 3, 4) × π/Ω(R). The bending mode decays rapidly due to horizontal mixing, while the breathing mode persists longer, driven by spiral arms. This difference of the two modes causes a transition to breathing dominance over time on the timescale of horizontal mixing (R/σR).

Fig. 7 First panel: spectrogram of the m = 1 bending mode. Lines indicate characteristic frequencies: Ω (solid), ± v (dashed), Ω ± v, and Ω ± 2 v (dotted). Second panel: spectrogram of the m = 2 breathing mode. Solid and dashed lines indicate Ω and Ω ± κ/2, respectively. Third panel: spectrogram of the surface m = 2 density wave. Solid and dashed lines indicate Ω and Ω ± κ/2, respectively. The bending mode shows distinct branches aligned with the resonant curves of Ω ± v and Ω−2 v. The breathing mode and the density wave rotate at nearly the same pattern speeds, suggesting that the breathing mode is excited by the bar and the tidally induced arms.

4.2 Spectral analysis

To investigate the connection between bending and breathing modes and density waves (such as bars and spiral arms), we performed a spectral analysis (Sellwood & Athanassoula 1986; Roškar et al. 2012; Chequers & Widrow 2017; Chequers et al. 2018; Poggio et al. 2021; Khachaturyants et al. 2022). This analysis provides insights into the frequencies of these modes across the galactocentric radius by visualising their power spectra in the radius-frequency space. We followed the method formalised by Chequers & Widrow (2017) for the computation of the spectrum.

The temporal discrete Fourier transform of the complex Fourier coefficients Q_n,m(R,t)(n ∈ bend, breath, dens) was defined as

$$F_{n,m}(R,{\omega }_k)=\sum _{j=0}^{N-1}{Q}_{n,m}(R,t_j)w(j)\exp \left(2\pi i\frac{jk}{N}\right),$$(14)

where w(j) is a window function. To suppress spectral leakage, we used the Gaussian function with a width of N/4. The time interval used is between the pericentre passage and the end of the simulation, with N = 185 time steps in this interval. The frequency ω_k is given by

$${\omega }_k=\frac{2\pi }{m}\frac{k}{N\mathrm{\Delta }t};k=-N/2\dots N/2,$$(15)

where Δ t = 9.78 Myr is the output cadence of the snapshot. The power spectrum was defined as

$$P_{n,m}(R,{\omega }_k)=\frac{1}{N\sum _{j}w(j)}{\left|F_{n,m}(R,{\omega }_k)\right|}^2.$$(16)

The first panel of Fig. 7 shows the power spectrum P_{n, m}(R, ω) on R – ω plane (spectrogram) of the m = 1 bending mode, which is the dominant Fourier component. Characteristic frequencies, such as Ω, ± v, Ω ± v, and Ω ± 2 v, are overplotted for comparison. At R ≳ 13 kpc, three distinct branches are visible: Ω ± v, corresponding to the pattern speed of the kinematic bending wave (Binney & Tremaine 2008; Chequers et al. 2018), and Ω−2 v, which has also been observed in other simulations (Chequers & Widrow 2017; Chequers et al. 2018; Poggio et al. 2021)⁴. Higher-order frequencies, such as Ω ± 3 v, are also present but significantly weaker. In the inner galaxy (R ∼ 5−13 kpc), a faint branch approximately follows Ω. Interestingly, this branch lies within the forbidden region (Ω−v < ω < Ω+v) according to WKB analysis (Binney & Tremaine 2008; Chequers et al. 2018). Its presence suggests that the linear approximation may break down in the inner galaxy, likely due to the stronger influence of self-gravity.

The second and third panels of Fig. 7 show the spectrograms of the m = 2 breathing and density components, respectively. Their similarities imply a strong connection between the breathing mode and non-axisymmetric density structures such as spiral arms and bars. A horizontal line at ω ∼ 40 km s⁻¹ kpc⁻¹ at R ≲ 4 kpc corresponds to the bar. Meanwhile, a branch along Ω−κ/2 is visible, aligning with the pattern speed of tidally induced spiral arms (Antoja et al. 2022). This matches the kinematic density wave pattern speed (Binney & Tremaine 2008).

In the outer galaxy (R ≳ 12 kpc), the spectrogram of the density mode shows a branch that follows Ω, consistent with dynamic spiral arms that corotate with disc stars (Baba et al. 2013). These dynamic arms differ from tidally induced spiral arms, which rotate more slowly relative to the disc. The correspondence between the m = 2 breathing mode and spiral arms supports the hypothesis that spiral structures play a significant role in driving the breathing mode.

Fig. 8 Maps of z–vz of the particles in the ‘solar neighbourhood’. The panels are colour-coded by number density (left column), density contrast (middle left column), radial velocity (middle right column), and azimuthal velocity (right column). The maps in the first row are at t = 0 Gyr. Second, third, and fourth rows correspond to 0.6 Gyr after the first pericentre passage, 0.5 Gyr after the second pericentre passage, and 0.5 Gyr after the third pericentre passage, respectively. Prominent one-arm phase spirals emerge in the solar neighbourhood after the first pericentre passage, becoming fainter after subsequent pericentre passages which have weaker impacts on the galactic disc.

5 Phase spirals

In the previous section, we focused on the global vertical oscillations of the disc, namely the bending and breathing modes, which provide insights into the large-scale impact of the dwarf galaxy. In this section, we shift our focus to the local scale, examining the formation and evolution of phase spirals-vertical phase-space features that are critical for understanding the disc’s response on smaller scales. These structures offer complementary insights into the dynamical processes studied in the global analysis.

5.1 Maps of z–v_z space in the solar neighbourhood

In Fig. 8, we show the z−v_z maps for the particles within 1 kpc from (x, y, z)=(−8.277, 0, 0) kpc, which corresponds to the solar position in the MW (GRAVITY Collaboration 2022). These maps highlight vertical phase-space structures and are divided into 90 × 90 bins, each with a size of 0.022 kpc × 1.33 km s⁻¹.

The left and middle-left panels are colour-coded by the particle count (n) and the density contrast (Δ ρ), respectively. The density contrast is defined as Δ ρ(z, v_z)=(n * G_1.5)(z, v_z) /(n *.G₃)(z, v_z), where n * G_σ represents the particle count convolved with the Gaussian filter with a standard deviation of σ pixels. This visualisation method, used in previous studies (e.g. Laporte et al. 2019; Antoja et al. 2023), enhances the visibility of phase spiral patterns. The middle-right and right panels are colour-coded by the deviations in mean radial velocity ($\mathrm{\Delta }{v}_R={\overline{v}}_R-{\overline{v}}_R\ast G_{10}$) and azimuthal velocity ($\mathrm{\Delta }{v}_{\varphi }=-{\overline{v}}_{\varphi }+{\overline{v}}_{\varphi }\ast G_{10}$) from the Gaussian-smoothed values, respectively.

The first row of Fig. 8 shows the z−v_z maps at t = 0 Gyr, before the dwarf’s pericentre passage. No asymmetric patterns are evident, indicating an initially symmetric phase-space distribution. While Bland-Hawthorn et al. (2019) identifies a quadrupole pattern in maps colour-coded by v_R due to the tilting of the velocity ellipsoid, this feature is removed here by subtracting the Gaussian-smoothed radial velocity.

In the second row, corresponding to t = 1.5 Gyr (0.6 Gyr after the first pericentre passage), prominent one-arm phase spirals are visible in maps colour-coded by Δ ρ, Δ v_R, and Δ v_ϕ. This one-arm feature is qualitatively consistent with the phase spiral observed in the solar neighbourhood (Antoja et al. 2018, 2023). The phase spiral in the maps colour-coded by radial and azimuthal velocities arises from the coupling between vertical and planar motion (Antoja et al. 2018; Binney & Schönrich 2018; Darling & Widrow 2019b; Bland-Hawthorn & Tepper-García 2021). The vertical orbital frequency depends not only on the vertical oscillation amplitude but also on the angular momentum (L_z=R v_ϕ). As a result, phase spirals are particularly pronounced in v_ϕ-coded maps, where angular momentum variations are highlighted. In contrast, v_R-coded spirals originate from the interplay of vertical and planar waves triggered by satellite perturbations (Antoja et al. 2023).

The third and fourth rows show the z−v_z maps at 0.5 Gyr after the second and third pericentre passages, respectively. Phase spirals remain visible, though they are notably fainter than those following the first pericentre passage. In the third row, spiral structures are more distinct in maps colour-coded by Δ v_R and Δ v_ϕ than in those colour-coded by Δ ρ. By contrast, the fourth row reveals faint spirals across all panels, suggesting significantly reduced perturbations. As shown in Fig. 1, the dwarf loses approximately 50% of its mass during each pericentre passage. This substantial mass loss weakens the perturbations, resulting in insufficient excitation of phase spirals during the third pericentre passage.

Fig. 9 Grid in the R−ϕ space. The grid size is 1 kpc × 20°. The contour map shows the normalised surface density at t = 1.47 Gyr.

5.2 Spatial variation of the phase spirals

In this section, we examine how the morphology of the phase spiral varies with the galactocentric radius and the azimuth. Fig. 9 shows the normalised surface density at t = 1.47 Gyr on the R−ϕ plane, where the galaxy rotates in the direction in which ϕ decreases (i.e. from top to bottom in the figure). We divided the R−ϕ space into a grid with a size of 1 kpc × 20° as shown in the figure and made the vertical phase-space maps for particles within 1 kpc of each grid point.

Fig. 10 shows the resulting phase-space maps, highlighting spatial variations in phase spirals at t = 1.47 Gyr. Each panel corresponds to a grid point in Fig. 9 and displays a vertical phase-space map colour-coded by the density contrast, Δ ρ. The horizontal and vertical axes are normalised by the standard deviations h_z and σ_z, respectively, and cover the range [−3.5, 3.5] ×[−3.5, 3.5]. To classify the phase spiral shapes, we performed a visual inspection: red panels indicate one-arm phase spirals, blue panels indicate two-arm phase spirals, and grey panels correspond to cases where the shapes are ambiguous.

Phase spirals or asymmetric distributions appear in most panels. In the outer galaxy, the structures are highly distorted, with some panels displaying U-shaped or C-shaped patterns rather than spirals. At R = 3 kpc, the spiral pattern is unclear due to the low density contrast. The phase spirals in the inner galaxy are more tightly wound, reflecting the faster vertical winding rate in this region. No clear trends in azimuthal variation are evident from Fig. 10. A more quantitative analysis-such as of parameterising features such as amplitude and rotation-is required to explore this further, as suggested by Alinder et al. (2023), but such an analysis is beyond the scope of this paper.

While one-arm phase spirals are widely observed, two-arm phase spirals are only found in specific regions, such as (R, ϕ)=(6 kpc,−30°). As discussed earlier, tidally induced spiral arms are a likely driver of the breathing mode, which may contribute to the formation of two-arm phase spirals⁵. This connection is particularly plausible in the inner galaxy, where the breathing mode dominates. However, Fig. 10 does not show a clear spatial correlation between the spiral arms and the two-arm phase spirals.

Fig. 11 shows the similar phase-space maps as Fig. 10, but for the snapshot at t = 1.65 Gyr. Compared to t = 1.47 Gyr (Fig. 10), two-arm phase spirals are more widely distributed, particularly between R = 6 kpc and 7 kpc. At R>8 kpc, the one-arm phase spirals are more widely observed than the two-arm phase spirals. This trend is consistent with the findings of Hunt et al. (2022), who reported that the two-arm phase spirals are observed only in the samples with small guiding radii. For further analysis, Appendix A.1 provides phase-space maps classified by guiding radius (R_g) and azimuthal angle (θ_ϕ) instead of R and ϕ.

Previous studies (Widrow et al. 2014; Banik et al. 2022, 2023) have shown that satellites can excite the breathing mode if their perturbations are impulsive, meaning that the satellites move sufficiently faster than disc stars. Analytical calculations by Banik et al. (2023) indicate that the bending mode dominates over the breathing mode in the solar neighbourhood and the inner galaxy; on the contrary, a substantial breathing mode amplitude is expected in the outer galaxy. Consistent with this prediction, in the N-body model by Hunt et al. (2021), two-arm phase spirals are observed only in the outer galaxy, and it is suggested that these phase spirals are associated with the breathing mode directly excited by external perturbations from the satellite. In our model, the two-arm phase spirals in the inner galaxy are likely excited by the spiral arms. In contrast, those in the outer galaxy, such as those observed around (R, ϕ)=(12 kpc,−50°) in Fig. 10 and (R, ϕ)=(11 kpc, 70°) in Fig. 11, might be directly excited by the satellite impact, consistent with the mechanism proposed in Hunt et al. (2021)’s model.

Some panels, such as (R, ϕ)=(6 kpc,−150°) in Fig. 10 and (9 kpc,−130°) in Fig. 11, exhibit ring-like structures rather than spirals. These features, absent in Gaia data, are of uncertain origin but may be linked to spiral arms, as they are located within arm regions.

Fig. 10 Variation of R−ϕ of the phase spiral in the snapshot at t= 1.47 Gyr. Each panel shows z/hz−vz/σz map for a set of particles grouped by R and ϕ. The maps are colour-coded by the density contrast, Δ ρ. Panels exhibiting one-arm phase spirals are painted in red tone, while those with two-arm phase spirals are painted in blue tone. Panels where the phase spiral shapes cannot be classified as either one-arm or two-arm are painted in grey tones. The surface density contour of Fig. 9 is overlaid. While one-arm phase spirals are widely observed, two-arm phase spirals are only found in limited areas.

Fig. 11 Same as Fig. 10 but for the snapshot at t = 1.65 Gyr. Two-arm phase spirals are more widely distributed than in Fig. 10, especially around R ∼ 6−7.5 kpc, consistent with observations showing two-arm spirals primarily for stars with small guiding radii.

5.3 Time evolution of the phase spiral

Finally, we examined the time evolution of the phase spirals. For all snapshots from t = 0.9 Gyr to 1.67 Gyr, we created z/h_z−v_z/σ_z maps at eighteen different positions at R = 8 kpc, equally spaced by 20° in azimuth. Fig. 12 illustrates these results, with rows representing time and columns representing azimuth. This type of visualisation was introduced by Bland-Hawthorn & Tepper-García (2021), who referred to it as a ‘chronogram’ in their paper. Each panel shows the phase-space map colour-coded by the density contrast, Δ ρ. The overlaid contours represent the normalised surface density, Σ(R = 8 kpc, ϕ, t)/Σ₀(R = 8 kpc, t), and are spaced by 0.1 between 1 and 1.8. Regions where Σ/Σ₀>1 are shaded to indicate arm regions.

At t ≲ 1.15 Gyr, phase spirals are not yet visible, but asymmetric density distributions are evident in the z/h_z−v_z/σ_z space. Around ϕ ∼−180°, phase spirals begin to emerge at t ∼ 1.15 Gyr (∼ 250 Myr after the first pericentre passage) and subsequently appear at other azimuths. By t ∼ 1.3 Gyr, phase spirals are visible throughout the entire ring at R = 8 kpc, though their prominence varies with azimuth. As expected, the phase spirals wind up over time, reflecting the underlying dynamics.

Two-arm phase spirals, highlighted in blue in Fig. 12, first appear at t ∼ 1.35 Gyr. These structures are limited to specific azimuthal ranges and shift positions following the galaxy’s rotation. They tend to occur more frequently in inter-arm regions. For example, a two-arm phase spiral emerges at ϕ = 30° at t ∼ 1.4 Gyr, moves leftward (in the direction of rotation), and becomes faint after entering the arm region at t ∼ 1.44 Gyr. Another example begins at ϕ=−10° at t ∼ 1.5 Gyr and disappears around ϕ=−100° at t ∼ 1.55 Gyr after entering an arm region. After entering the arm region, these two-arm spirals transform into ring-like structures, as seen at ϕ ∼−150° at t ∼ 1.57 Gyr. This transformation suggests that the breathing perturbation from spiral arms deforms phase spirals into ring structures.

Fig. 13 shows the time evolution of phase spirals at R= 6 kpc. Phase spirals appear at t ∼ 1.15 Gyr and become increasingly tightly wound, consistent with observations at R = 8 kpc. However, two-arm phase spirals emerge earlier at t ∼ 1.25 Gyr, reflecting the faster transition from bending mode to breathing mode dominance in the inner galaxy. Unlike at R = 8 kpc, twoarm phase spirals at R = 6 kpc are not confined to specific azimuths. At earlier times (t ≲ 1.5 Gyr), the patterns resemble those at R = 8 kpc. After t ∼ 1.5 Gyr, particularly around t ∼ 1.65 Gyr, two-arm phase spirals become more widespread, likely due to weakening spiral arms. Between t = 1.0 Gyr and 1.6 Gyr, the spiral arm density at R = 6 kpc decreases by ∼ 50%, blurring the boundaries between arm and inter-arm regions. Additionally, the bar may influence phase spirals at R = 6 kpc, as this radius is near the bar’s corotation radius.

At both R = 6 kpc and R = 8 kpc, two-arm phase spirals appear intermittently rather than continuously after their initial emergence. This temporal irregularity is likely due to the oscillation in spiral strength and the breathing amplitude (see Fig. 4). Although predicting the exact timing and location of two-arm spirals is challenging, their first appearance correlates with the transition to breathing mode dominance.

Comparing Figs. 12, 13, and the third panel of Fig. 5, we find that the two-arm phase spiral emerges ∼ 200 Myr after the dominant mode transition at both R = 6 kpc and 8 kpc. This trend holds for other radii between R = 5 kpc and R = 10 kpc, where two-arm spirals consistently appear 200–250 Myr after the mode transition. The delay in the emergence of the two-arm phase spiral relative to the dominant mode transition occurs because the bending and breathing amplitudes reflect the current states of these oscillation modes, while the phase spirals correspond to their histories.

In Appendix A.2, we present the time evolution of the vertical phase space maps colour-coded by Δ v_R and Δ v_ϕ. The phase spirals labelled by Δ v_R and Δ v_ϕ emerge earlier than those labelled by Δ ρ, and their morphologies (one-arm or two-arm) do not always coincide.

Fig. 12 Time evolution of the phase spiral at R = 8 kpc. Columns represent azimuths equally spaced by 20°. Rows correspond to time, from t = 0.9 Gyr to 1.28 Gyr on the left block, and t = 1.29 Gyr to 1.67 Gyr on the right block. Each panel shows z/hz−vz/σz maps colour-coded by the density contrast, Δ ρ. Panels showing two-arm phase spirals are painted in blue tone. Contours indicate the normalised surface density, Σ(R = 8 kpc, ϕ, t)/Σ0(R = 8 kpc, t). Arm regions (Σ/Σ0>1) are shaded. One-arm spirals emerge first, followed by intermittent appearances of two-arm spirals around 200-250 Myr after the bending-to-breathing mode transition.

Fig. 13 Same as Fig. 12 but at R = 6 kpc. Two-arm phase spirals emerge earlier at R = 6 kpc than at R = 8 kpc, reflecting the faster transition to breathing mode dominance in the inner galaxy.

6 Discussions

6.1 Dating the perturbation to the MW disc

The phase spiral contains information about the timing of the perturbation event, as it winds up over time at a fixed rate determined by the galactic potential. Several previous studies (e.g. Antoja et al. 2018, 2023; Frankel et al. 2023; Darragh-Ford et al. 2023) estimated the dynamical age of the one-arm phase spiral to be ∼ 200−900 Myr by ‘unwinding’ the spiral under assumed galactic potentials. However, these analyses often neglect the disc’s self-gravity, which slows the winding rate of phase spirals (Widrow 2023) and could lead to an underestimation of their dynamical age.

An alternative approach was proposed by Antoja et al. (2022), who linked the perturbation time to tidally induced spiral arms. In our simulation, tidal forces from a satellite induce spiral arms in the disc that rotate at Ω−κ/2. These arms create a global velocity pattern that manifests as a wave in v_R as a function of L_z when observed within a narrow azimuthal range. The wavelength of this wave corresponds to the time since the formation of the arms, as the wave shortens over time due to the winding of the arms. Friske & Schönrich (2019) detected a similar L_z−v_R wave in the Gaia data, and Antoja et al. (2022) identified two distinct wavelengths corresponding to the arm formation times of < 0.6 Gyr and 0.8-2.1 Gyr through Fourier analysis. Although the origin of these waves remains debated (Cao et al. 2024), this method offers a complementary avenue for estimating perturbation times.

Our results suggest a new method for estimating the perturbation time using the relative strengths of the bending and breathing modes. In our simulations, we observed a transition from bending-dominated to breathing-dominated vertical oscillations. This transition progresses outward from the inner to the outer galaxy, with the timescale determined by R /σ_R. While it is challenging to determine azimuthally averaged amplitude ratios from observational data as in Section 3.2, local observations of phase spirals can provide constraints. For example, in Section 5.3, we showed that two-arm phase spirals emerge ∼ 200−250 Myr after the bending-to-breathing transition. The detection of two-arm phase spirals in Gaia data implies that part of the MW disc transitioned to a breathing-dominated state at least ∼ 200 Myr ago.

The two-arm phase spiral is predominantly observed for stars with small guiding radii, roughly R_g ≲ 7 kpc (Hunt et al. 2022; Alinder et al. 2024). The bending-to-breathing transition takes R/σ_R ∼ 200 Myr at R = 7 kpc, where σ_R ∼ 35 km s⁻¹ (Bland-Hawthorn & Gerhard 2016; Gaia Collaboration 2018b). This suggests that the MW disc was perturbed more than ∼ 400 Myr ago (200 Myr for the transition plus 200 Myr for the emergence of two-arm phase spirals). This estimate is consistent with the previous studies (Antoja et al. 2018, 2022, 2023; Li & Widrow 2021, 2023; Frankel et al. 2023; Darragh-Ford et al. 2023).

6.2 Sgr mass problem

Previous studies (e.g. Binney & Schönrich 2018; Bennett et al. 2022) suggest that satellites must be heavier than ∼ 3 × 10¹⁰ M_⊙ to excite phase spirals. This required mass is several orders of magnitude greater than the present-day mass of Sgr, estimated at ∼ 5 × 10⁸ M_⊙ (Vasiliev & Belokurov 2020). In our dwarf model, the satellite initially has a total mass of 5 × 10¹⁰ M_⊙ and loses 50% of its mass during each orbital loop, reducing its final mass to match the present-day Sgr. However, the phase spirals produced in the final epoch of the simulation are much fainter than those observed in the Gaia data.

This discrepancy between the theoretically required mass and the observationally estimated mass may be resolved by properly accounting for Sgr’s mass loss history. Bland-Hawthorn & Tepper-García (2021) suggest that the Sgr could have excited the phase spiral during earlier pericentre passages if it had experienced mass loss of 0.5−1 dex per orbit loop. Such high mass loss rates are dynamically plausible. Jiang & Binney (2000), using a semi-analytic model and N-body simulations, demonstrated that a wide range of orbital histories can reproduce the present-day position and velocity of the Sgr. An extreme scenario involves Sgr initially having a mass of ∼ 10¹¹ M_⊙ and being located at r ∼ 250 kpc from the MW centre. By fine-tuning the initial conditions of the Sgr model, it might be possible to simultaneously reproduce the phase spiral and the present-day state of Sgr. Such adjustments could help resolve the tension between the mass required to excite the phase spiral and Sgr’s currently observed mass.

6.3 Comparison with an isolated galaxy model

Breathing modes are not exclusive to perturbed galaxies; isolated galaxies can also exhibit breathing modes if they possess spiral arms or bars. Asano et al. (2024) investigate the relation between the breathing mode and the evolution of dynamic spiral arms in an isolated galaxy model. In this section, we compare the characteristics of breathing modes in an isolated galaxy with those in a perturbed galaxy, focusing on their spatial alignment and amplitude.

The upper panel of Fig. 14 shows the breathing velocity, v_{z, > 0}−v_{z, < 0}, around a spiral arm accompanying strong breathing mode in the isolated model of Fujii et al. (2019). In this frame, the galaxy rotates from top to bottom. The contours represent the normalised surface density, Σ(R, ϕ)/Σ₀(R), and the black dots indicate its peaks, detected using the signal.find_peaks routine in scipy package (Virtanen et al. 2020). The compressing breathing mode (indicated by the blue colour) is observed within the spiral arms, with no systematic displacement between the compressing area and the density peaks.

The lower panel of Fig. 14 shows the same plot for the perturbed model. In this case, the breathing signature is observed globally along the tidally induced spiral arm. The compressing mode is stronger on the trailing side of the arm than on the leading side, and the expanding mode is more prominent in the inter-arm regions.

This displacement is more evident in Fig. 15, where the breathing velocity and the normalised surface density are shown as functions of ϕ, averaged over R between R = 8 kpc and R = 8.5 kpc. In the isolated model (upper panel), the minimum of the breathing velocity almost coincides with the maximum of the density at ϕ ∼ 190°. In contrast, in the perturbed model (lower panel), the density peaks are located at ϕ ∼ 20°, while the minimum breathing velocity is shifted to ϕ ∼ 35°.

The difference in the spatial alignment of the breathing modes between the two models can be attributed to the distinct pattern speeds of dynamic arms versus tidally induced arms. For density-wave-like arms with pattern speeds slower than the galaxy’s rotation, compressing breathing modes and expanding breathing modes appear on the trailing and leading sides of the arms, respectively (Debattista 2014; Faure et al. 2014; Monari et al. 2016a). This scenario aligns with the tidally induced arms in our perturbed model, where the breathing velocity minimum is shifted to the trailing side. Conversely, dynamic arms, which corotate with disc stars at every radius, show no displacement between arm locations and breathing modes (Baba et al. 2013).

The amplitude of the breathing mode also differs between the models. In the perturbed model, the amplitude is higher than in the isolated model, as tidally induced arms are stronger and longer-lived than dynamic arms. This difference in amplitude affects the morphology of the phase spiral. Fig. 16 compares the two-arm phase spirals in the isolated and perturbed models, showing that the phase spiral in the isolated model is fainter.

Asano et al. (2024) conclude that the Local arm aligns with the dynamic arm model because there is no systematic displacement between the compressing breathing mode and the arm location. They also discussed the growth and disruption of the other arms based on the tentative signals of the breathing mode around them, assuming the dynamic arm model. However, they do not rule out the tidally induced arm model, as the spatial correlations between the arms and the breathing mode are less clear compared to the Local arm due to the observational uncertainties. The detection of the prominent two-arm phase spiral suggests that the MW has strong spiral arms, even if not tidally induced arms.

Fig. 14 Breathing signature around spiral arms. Upper panel: breathing signature in the isolated model. Colour map shows the breathing velocity, vz, > 0−vz, > 0. Contours indicate the normalised density, Σ(R, ϕ)/Σ0(R). Black dots indicate the density peaks. Lower panel: same as the upper panel but in the perturbed model. In the perturbed model (tidally induced arms), the breathing velocity minimum is displaced from the density peak to the trailing side, unlike in the isolated model (dynamic arms) where they are aligned.

Fig. 15 Breathing velocity and normalised density as functions of the azimuth. Upper panel: isolated model. The solid line and dashed line represent the breathing velocity and the normalised density, respectively, averaged over R between R = 8 kpc and 8.5 kpc. The arrow indicates the direction of the galaxy rotation. Lower panel: same as the upper panel but in the perturbed model. This plot quantifies the displacement in Fig. 14, showing that the breathing velocity minimum is shifted to the trailing side relative to the density peak in the perturbed model with tidally induced arms. This contrasts with the alignment seen in the isolated model with dynamic arms.

Fig. 16 Two-arm phase spiral in the isolated model (left) and in the perturbed model (right). The two-arm phase spiral is fainter in the isolated model than in the perturbed model.

6.4 Model limitations

Our high-resolution cold disc model captures the spatial and temporal evolution of vertical oscillations in the MW disc in detail. However, several important components that could influence the disc’s dynamical response are not included. While the MW disc is composed of multiple components, specifically the thin disc, thick disc, and gas disc, our model only has a single component corresponding to the thin disc. Observational estimates suggest that the thick and gas discs each account for only ∼ 10−20% of the total disc mass (Bland-Hawthorn & Gerhard 2016; Nakanishi & Sofue 2016; McMillan 2017). Since the thin disc makes the dominant contribution to the total disc potential, it is reasonable to assume that the thin disc dominates the overall dynamical response to external perturbations. Nonetheless, we discuss below the potential influence of the thick and gas discs based on insights from recent studies.

Stelea et al. (2024) performed N-body simulations of warm and cold stellar discs perturbed by the Sgr and the Large Magellanic Cloud (LMC). Their models exhibit only minor differences in the resulting bending amplitudes, implying that the inclusion of a thick disc (i.e. a dynamically hotter component) does not significantly alter the global bending motion. However, the two models show notable difference in the v_R fields, which are related to the tidally induced arms. This suggests that the thick disc may modulate the amplitude of breathing modes, which are excited by the spiral arms.

The impact of the gas disc is less well understood. Tepper-García et al. (2022) study the evolution of vertical corrugations in a combined stellar+gas disc using hybrid N-body and hydrodynamical simulations. They demonstrate that, although stellar and gas corrugations are initially excited in phase, the gas component is damped more quickly due to dissipative processes. Based on their results, if the interaction between the gas disc and the stellar disc is not negligible, gas disc may influence the evolution of bending modes in the stellar component. We also expect that self-gravity is more important in stellar+gas discs than in pure stellar discs, as the gas helps keep the disc dynamically cold. Therefore, high-resolution N-body and hydrodynamical simulations are important to understand influences of self-gravity on the disc oscillation. Additionally, if we can include star formations in such simulations, they will allow us to test the possible correlation between the Sgr impact and the star formation activity in the MW disc suggested by Ruiz-Lara et al. (2020).

Our study focuses on the perturbation of the Sgr on the vertical structure of the MW disc. However, it is important to recognise that the MW is also interacting with the LMC, whose gravitational influence cannot be neglected. Some recent studies (e.g. Laporte et al. 2018; Stelea et al. 2024) performed N-body simulations of the MW-Sgr-LMC system to investigate the combined impact of these two satellites. One of the notable features seen in their models is the disc warp similar to the warp in the MW disc. Although our Sgr-only model also produces a warp, its orientation in the present-day snapshot (t = 1.78 Gyr) does not align with that of the observed MW. In our model, the outer disc exhibits a vertical displacement of z < 0 at y > 0 and z > 0 at y < 0, whereas the observed MW warp displays the opposite pattern. This discrepancy suggests that the LMC plays a crucial role in shaping the actual Galactic warp.

The outer disc appears to be strongly perturbed by the combined influences of Sgr and the LMC, while the inner disc is thought to be primarily affected by Sgr. The LMC is currently located at a Galactocentric distance of ∼ 50 kpc, having recently passed its pericentre. The disc is only beginning to respond to the LMC’s perturbation from the outer part. There has not yet been sufficient time for vertical features such as corrugations or phase spirals to fully develop (Laporte et al. 2019). Therefore, our Sgr-only model remains valuable for studying the disc dynamics at inner and intermediate radii, where the Sgr’s influence dominates. As an example, we used the ratio R/σ_R to trace the spatial transition from bending-dominated to breathingdominated regimes. This diagnostic should remain valid for R ≲ 15 kpc, even in the presence of the LMC. At larger radii, however, models incorporating both Sgr and the LMC are necessary for a more accurate representation.

In addition to its direct perturbation, the LMC may also influence the orbit of the Sgr dwarf, thereby altering its interaction history with the MW. Vasiliev et al. (2021) demonstrate that including the LMC in N-body simulations of the MW-Sgr interaction modifies Sgr’s orbits, particularly its eccentricity. This affects the pericentre distances and potentially alters the strength of Sgr’s impact during previous pericentre passages.

In summary, Sgr appears to be the primary perturber of the MW disc, at least in the solar neighbourhood, and the response of the thin disc plays a key role in determining the overall dynamical behaviour. However, a comprehensive understanding of the MW’s response to external perturbations requires simulations of multi-component discs that include both Sgr and the LMC, and we plan to perform such simulations in future work. They will help us reconstruct the history of perturbation events encoded in the bending and breathing signatures in the MW disc.

7 Conclusions

In this study, we used high-resolution N-body simulations to investigate the interaction between the MW and the Sgr dwarf galaxy. This interaction generates vertical oscillations in the MW disc, including bending and breathing modes, as well as phase spirals associated with these oscillations. Our findings highlight the interplay between direct satellite perturbations and indirect effects mediated by tidally induced spiral arms, shaping the vertical structure of the MW disc. Our main results are summarised as follows:

1. The satellite’s perturbation directly excites the bending mode across a wide radius in the galactic disc. Simultaneously, it induces the spiral arms within the disc;

2. The tidally induced spiral arms subsequently excite the breathing mode in the disc. As such, the breathing mode is the indirect result of the satellite interaction;

3. Initially, the bending mode dominates, but the dominant mode transitions to the breathing mode over time. This transition originates from the different timescales of the two modes: the bending mode decays rapidly mainly due to horizontal mixing, while the breathing mode persists, sustained by the spiral arms. The transition occurs faster in the inner galaxy than in the outer regions;

4. The bending mode is dominated by its m = 1 Fourier component, with its spectrogram displaying distinct branches aligned with the resonant curves of Ω ± v and Ω−2 v. In contrast, the dominant Fourier component for the breathing mode is m = 2, which rotates at ∼ 40, km s⁻¹ kpc⁻¹ in the inner galaxy and Ω−κ/2 in the outer galaxy, corresponding to the pattern speeds of the bar and the tidally induced spiral arms, respectively;

5. The simulation successfully reproduces the one-arm phase spiral observed in the solar neighbourhood. Additionally, it reveals two-arm phase spirals, particularly in the inner galaxy, linked to the breathing mode excited by spiral arms. The simulations also indicate that the two-arm phase spiral emerges approximately 200–250 Myr after the transition to the breathing-dominated state;

6. The detection of the two-arm phase spiral in the Gaia data suggests that parts of the MW disc transitioned to a breathing-dominated state at least ∼ 200 Myr ago. Considering the bending-to-breathing transition timescale, we estimate that the MW disc was perturbed more than ∼ 400 Myr ago.

These results provide insights into the MW disc’s response to external perturbations, particularly from the Sgr dwarf galaxy, and highlight the importance of both direct and indirect effects in shaping the disc’s vertical structure. We plan to extend this work by incorporating the LMC and a multi-component disc model to further investigate the MW’s response to external perturbations and to reconstruct the history of perturbation events encoded in the MW disc.

Acknowledgements

We thank the anonymous referee for the useful comments. We appreciate discussions with the members of the Gaia-UB group. This research used computational resources of Pegasus and Cygnus provided by Multidisciplinary Cooperative Research Program in Center for Computational Sciences, University of Tsukuba. Part of the numerical computations were carried out on HA-PACS at Tsukuba University, Piz Daint at the Swiss National Supercomputing Centre and Little Green Machine II. In this research work, we used the “mdx: a platform for the data-driven future”. TA was supported by Japan Society for the Promotion of Science (JSPS) Research Fellowships for Young Scientists and JSPS Overseas Challenge Program for Young Researchers. This work was supported by Grant-in-Aid for JSPS Fellows Number JP22J11943. JB acknowledges support from JSPS under Grant Numbers 21 K 03633 and 21H00054. We acknowledge the grants PID2021-125451NA-I00 and CNS2022-135232 funded by MICIU/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, by the “European Union” and by the “European Union Next Generation EU/PRTR”. This work made use of the following software packages: Bonsai (Bédorf et al. 2012, 2014), Galactics (Kuijken & Dubinski 1995; Widrow & Dubinski 2005; Widrow et al. 2008), astropy (Astropy Collaboration 2013, 2018, 2022), Jupyter (Perez & Granger 2007; Kluyver et al. 2016), matplotlib (Hunter 2007), numpy (Harris et al. 2020), pandas (Wes McKinney 2010; The pandas development team 2024), scipy (Virtanen et al. 2020; Gommers et al. 2023) and Agama (Vasiliev 2019). This research has made use of NASA’s Astrophysics Data System. Software citation information aggregated using The Software Citation Station (Wagg & Broekgaarden 2024; Wagg et al. 2024).

Appendix A Additional material about the phase spiral

A.1 Dependence of R_g and θ_ϕ of the phase spiral

We created vertical phase-space maps by dividing the local sample based on the guiding radius (R_g) and azimuthal angle (θ_ϕ) in a similar way to how Hunt et al. (2021, 2022) did. This sample division allowed us to estimate the global phase mixing state from the local sample (Li 2021; Gandhi et al. 2022). We first constructed the potential model from the N-body snapshot using Agama (Vasiliev 2019). We applied multipole expansion to the DM halo and the classical bulge and azimuthal harmonic expansion to the disc. We converted the angular momentum, L_z, to the guiding radius, R_g, by interpolating the R–L_z relation obtained from the potential model. The azimuthal angle was computed with Stäckel fudge (Binney 2012) implemented in Agama. It is important to note that the approximation of axisymmetric potentials may lead to significant biases in the estimation of angle-action variables when the galaxies have strong spiral arms (Debattista et al. 2025).

Fig. A.1 shows the face-on map of the galactic disc at t= 1.65 Gyr, where six circles indicate the local regions labelled A-F. Each circle is centred at R = 9 kpc and has a radius of 1.5 kpc. Fig. A.2 shows the R_g−θ_ϕ dependence of the phase spiral in Region D. Each panel shows the z/h_z−v_z/σ_z map colour-coded by Δ ρ. Two-arm phase spirals are observed in the small R_g (inner galaxy) bins, while one-arm phase spirals appear in the large R_g (outer galaxy) bins. This spatial variation in spiral morphology consistent with what is observed in the Gaia data (Hunt et al. 2022). Figs. A.3–A. 7 show the same plots for the other regions. Whereas Region A (Fig. A.3) shows a similar trend, the twoarm and one-arm phase spirals are not clearly separated by R_g in the Regions B, C, E and F. In these regions, the two-arm phase spirals tend to be observed in the intermediate R_g range around ∼ 7−9 kpc rather than the inner galaxy.

Fig. A.1 Face-on density map at t = 1.65 Gyr. Six circles indicate regions A−F.

Fig. A.2 Dependence of Rg−θϕ of the phase spiral in Region D. Each panel shows the z/hz−vz/σz map colour-coded by Δ ρ. The galaxy rotates from top to bottom in this frame.

Fig. A.3 Same as Fig. A.2 but for Region A.

Fig. A.4 Same as Fig. A.2 but for Region B.

Fig. A.5 Same as Fig. A.2 but for Region C.

Fig. A.6 Same as Fig. A.2 but for Region E.

Fig. A.7 Same as Fig. A.2 but for Region F.

Fig. A.8 Time evolution of the phase spiral at R = 8 kpc. The z/hz−vz/σz maps are shown in the same way as Fig. 12 but colour-coded by Δ vR instead of Δ ρ.

A.2 Time evolution of phase spirals label by v_R and v_ϕ

Fig. A.8 and Fig. A.9 show the time evolution of the phase spiral at R = 8 kpc. They are visualised similarly to Fig. 12 but are colour-coded by Δ v_R and Δ v_ϕ, respectively. One-arm phase-spirals start to emerge around t ∼ 1.05 Gyr in the maps colour-coded by Δ v_R (Fig. A.8) and Δ v_ϕ (Fig. A.9), whereas they emerge around t ∼ 1.15 Gyr in the maps colour-coded by Δ ρ (Fig. 12). Morphologies (i.e. one-arm or two-arm) of the phase spirals labelled by Δ ρ, Δ v_R, and Δ v_ϕ do not always coincide, suggesting a complex vertical phase mixing process coupled with planar motion.

Fig. A.9 Same as Fig. A.8 but colour-coded by Δ vϕ.

References

All Tables

All Figures

Fig. 1 Upper panel: position of the dwarf as a function of time. Blue solid, orange dashed, and green dotted lines show $r=\sqrt{x^2+y^2+z^2},R=\sqrt{x^2+y^2}$, and z, respectively. Lower panel: mass of the dwarf as a function of time. Blue solid and orange dashed lines show the total mass and stellar mass enclosed within 5 kpc. The Green dotted line shows the stellar mass fraction. The red vertical line indicates the time of t = 1.78 Gyr when the dwarf’s position is close to the present-day position of the Sgr. The final total mass is consistent with the presentday Sgr mass.

In the text

Fig. 2 Projection of the dwarf’s orbit in the x–y plane (upper panel) and in the x–z plane (lower panel). The cyan line indicates the orbital trajectory of the dwarf. The cyan point and cyan solid arrow indicate the dwarf’s position and velocity at t = 1.78 Gyr, respectively. The red square and red dashed arrow indicate those for the Sgr at the present day estimated by Vasiliev & Belokurov (2020). Purple dots show the distribution of the stellar particles stripped away from the dwarf. The surface density of the host disc and the solar position are displayed on the background. The dwarf follows a nearly polar orbit similar to the Sgr, leaving a stellar stream of stripped particles along its orbital path.

In the text

	Fig. 3 Face-on maps of the disc at t = 1.15 Gyr, which is ∼ 250 Myr after the first pericentre passage of the dwarf. Tidal forces from the dwarf induce prominent two-arm spiral arms and cause vertical corrugations (bending) and a breathing pattern in the galactic disc. The breathing pattern is aligned with the spiral arms.
In the text

Fig. 4 Time evolution of the bending mode, breathing mode, and surface density in the ring of R = 8–8.5 kpc. First panel: total amplitudes of the bending (solid line) and breathing (dashed line) as functions of time. Second panel: Fourier amplitudes of the bending mode for m = 0, 1, 2, 3, and 4. Different line styles correspond to different Fourier components as indicated in the legend. Third panel: same as the second panel but for the breathing mode. Fourth panel: time evolution of the Fourier amplitudes of the surface density. In all panels, vertical lines indicate the timing of the pericentre passages of the dwarf. The bending mode responds sharply to pericentre passages and decays quickly, while the breathing mode increases more slowly and persists longer, leading to a transition from bending-dominated to breathing-dominated state. The m = 1 and m = 2 are the dominant Fourier components for the bending and breathing modes, respectively.

In the text

Fig. 5 Bending amplitude, breathing amplitude, and their ratio as functions of time and galactocentric radius. First panel: total amplitude of the bending mode. Second panel: total amplitude of the breathing mode. Third panel: amplitude ratio between the bending mode and the breathing mode. The vertical lines indicate the times of the dwarf’s pericentre passages. The dominant mode changes from bending to breathing. This transition occurs faster in the inner galaxy.

In the text

Fig. 6 Evolution of the bending mode and the breathing mode in the single-impact model. First panel: bending and breathing amplitude ratio as a faction of time and the galactocentric radius. The horizontal axis indicates the time since the pericentre passage of the dwarf. A vertical dotted line indicates the pericentre time t–t0 = 0. Dashed curve represents R/σR(R). Second panel: bending amplitude as a function of time and the galactocentric radius. The amplitude is normalised by the maximum amplitude at each radius. Solid curves represent (1, 3, 5, 7) × π/v(R)−50 Myr. Third panel: same as the second panel but for the breathing mode. Solid lines represent (1, 2, 3, 4) × π/Ω(R). The bending mode decays rapidly due to horizontal mixing, while the breathing mode persists longer, driven by spiral arms. This difference of the two modes causes a transition to breathing dominance over time on the timescale of horizontal mixing (R/σR).

In the text

Fig. 7 First panel: spectrogram of the m = 1 bending mode. Lines indicate characteristic frequencies: Ω (solid), ± v (dashed), Ω ± v, and Ω ± 2 v (dotted). Second panel: spectrogram of the m = 2 breathing mode. Solid and dashed lines indicate Ω and Ω ± κ/2, respectively. Third panel: spectrogram of the surface m = 2 density wave. Solid and dashed lines indicate Ω and Ω ± κ/2, respectively. The bending mode shows distinct branches aligned with the resonant curves of Ω ± v and Ω−2 v. The breathing mode and the density wave rotate at nearly the same pattern speeds, suggesting that the breathing mode is excited by the bar and the tidally induced arms.

In the text

Fig. 8 Maps of z–vz of the particles in the ‘solar neighbourhood’. The panels are colour-coded by number density (left column), density contrast (middle left column), radial velocity (middle right column), and azimuthal velocity (right column). The maps in the first row are at t = 0 Gyr. Second, third, and fourth rows correspond to 0.6 Gyr after the first pericentre passage, 0.5 Gyr after the second pericentre passage, and 0.5 Gyr after the third pericentre passage, respectively. Prominent one-arm phase spirals emerge in the solar neighbourhood after the first pericentre passage, becoming fainter after subsequent pericentre passages which have weaker impacts on the galactic disc.

In the text

	Fig. 9 Grid in the R−ϕ space. The grid size is 1 kpc × 20°. The contour map shows the normalised surface density at t = 1.47 Gyr.
In the text

Fig. 10 Variation of R−ϕ of the phase spiral in the snapshot at t= 1.47 Gyr. Each panel shows z/hz−vz/σz map for a set of particles grouped by R and ϕ. The maps are colour-coded by the density contrast, Δ ρ. Panels exhibiting one-arm phase spirals are painted in red tone, while those with two-arm phase spirals are painted in blue tone. Panels where the phase spiral shapes cannot be classified as either one-arm or two-arm are painted in grey tones. The surface density contour of Fig. 9 is overlaid. While one-arm phase spirals are widely observed, two-arm phase spirals are only found in limited areas.

In the text

	Fig. 11 Same as Fig. 10 but for the snapshot at t = 1.65 Gyr. Two-arm phase spirals are more widely distributed than in Fig. 10, especially around R ∼ 6−7.5 kpc, consistent with observations showing two-arm spirals primarily for stars with small guiding radii.
In the text

Fig. 12 Time evolution of the phase spiral at R = 8 kpc. Columns represent azimuths equally spaced by 20°. Rows correspond to time, from t = 0.9 Gyr to 1.28 Gyr on the left block, and t = 1.29 Gyr to 1.67 Gyr on the right block. Each panel shows z/hz−vz/σz maps colour-coded by the density contrast, Δ ρ. Panels showing two-arm phase spirals are painted in blue tone. Contours indicate the normalised surface density, Σ(R = 8 kpc, ϕ, t)/Σ0(R = 8 kpc, t). Arm regions (Σ/Σ0>1) are shaded. One-arm spirals emerge first, followed by intermittent appearances of two-arm spirals around 200-250 Myr after the bending-to-breathing mode transition.

In the text

	Fig. 13 Same as Fig. 12 but at R = 6 kpc. Two-arm phase spirals emerge earlier at R = 6 kpc than at R = 8 kpc, reflecting the faster transition to breathing mode dominance in the inner galaxy.
In the text

Fig. 14 Breathing signature around spiral arms. Upper panel: breathing signature in the isolated model. Colour map shows the breathing velocity, vz, > 0−vz, > 0. Contours indicate the normalised density, Σ(R, ϕ)/Σ0(R). Black dots indicate the density peaks. Lower panel: same as the upper panel but in the perturbed model. In the perturbed model (tidally induced arms), the breathing velocity minimum is displaced from the density peak to the trailing side, unlike in the isolated model (dynamic arms) where they are aligned.

In the text

Fig. 15 Breathing velocity and normalised density as functions of the azimuth. Upper panel: isolated model. The solid line and dashed line represent the breathing velocity and the normalised density, respectively, averaged over R between R = 8 kpc and 8.5 kpc. The arrow indicates the direction of the galaxy rotation. Lower panel: same as the upper panel but in the perturbed model. This plot quantifies the displacement in Fig. 14, showing that the breathing velocity minimum is shifted to the trailing side relative to the density peak in the perturbed model with tidally induced arms. This contrasts with the alignment seen in the isolated model with dynamic arms.

In the text

	Fig. 16 Two-arm phase spiral in the isolated model (left) and in the perturbed model (right). The two-arm phase spiral is fainter in the isolated model than in the perturbed model.
In the text

	Fig. A.1 Face-on density map at t = 1.65 Gyr. Six circles indicate regions A−F.
In the text

	Fig. A.2 Dependence of Rg−θϕ of the phase spiral in Region D. Each panel shows the z/hz−vz/σz map colour-coded by Δ ρ. The galaxy rotates from top to bottom in this frame.
In the text

	Fig. A.3 Same as Fig. A.2 but for Region A.
In the text

	Fig. A.4 Same as Fig. A.2 but for Region B.
In the text

	Fig. A.5 Same as Fig. A.2 but for Region C.
In the text

	Fig. A.6 Same as Fig. A.2 but for Region E.
In the text

	Fig. A.7 Same as Fig. A.2 but for Region F.
In the text

	Fig. A.8 Time evolution of the phase spiral at R = 8 kpc. The z/hz−vz/σz maps are shown in the same way as Fig. 12 but colour-coded by Δ vR instead of Δ ρ.
In the text

	Fig. A.9 Same as Fig. A.8 but colour-coded by Δ vϕ.
In the text