<figure><img src="image_1.png" /></figure>

# A Potential Dynamical Origin of the Galactic Disk Warp: The Gaia–Sausage–Enceladus Major Merger

Mingji Deng<fnref n="1" />, Cuihua Du<fnref n="1" />, Yanbin Yang<fnref n="2" />, Jiwei Liao<fnref n="1" />, and Dashuang Ye<fnref n="1" />

Received 2024 April 8; revised 2024 August 30; accepted 2024 September 3; published 2024 October 23

## Abstract

Previous studies have revealed that the Galactic warp is a long-lived, nonsteady, and asymmetric structure. There is a need for a model that accounts for the warp's long-term evolution. Given that this structure has persisted for over 5 Gyr, its timeline may coincide with the completion of the Gaia–Sausage–Enceladus (GSE) merger. Recent studies indicate that the GSE, the significant merger of our Galaxy, was likely a gas-rich merger and the large amount of gas introduced could have created a profound impact on the Galactic morphology. This study utilizes GIZMO simulation code to construct a gas-rich GSE merger. By reconstructing the observed characteristics of the GSE, we successfully reproduce the disk warp and capture nearly all of its documented features, which align closely with observational data from both stellar and gas disks. This simulation demonstrates the possibility that a single major merger could generate the Galactic warp amplitude and precession. Furthermore, the analysis of the warp's long-term evolution may offer more clues into the formation history of the Milky Way.

Unified Astronomy Thesaurus concepts: Galaxy structure (622); Milky Way disk (1050); Milky Way dark matter halo (1049); Galaxy mergers (608); Hydrodynamical simulations (767)

## 1. Introduction

The most common asymmetrical structure in many disk galaxies is the disk warp, as documented in previous studies (M. L. Sánchez-Saavedra et al. 1990, 2003; V. Reshetnikov & F. Combes 1998). The Milky Way (MW), a typical example of disk and spiral galaxies, also exhibits a clear disk warp, confirmed by various works (F. J. Kerr 1957; H. T. Freudenreich et al. 1994; E. S. Levine et al. 2006; X. Chen et al. 2019). Many mechanisms of the Galactic warp have been proposed, including the inflow of intergalactic matter into the halo (E. C. Ostriker & J. J. Binney 1989; T. Quinn & J. Binney 1992; I.-G. Jiang & J. Binney 1999), inflow directly onto the Galactic disk (M. López-Corredoira et al. 2002), magnetic fields that exist between galaxies (E. Battaner & J. Jiménez-Vicente 1998), the interaction between satellite galaxies like Sagittarius (J. Bailin 2003) or Magellanic Clouds (M. D. Weinberg & L. Blitz 2006) with the disk, and the bending instability and self-excited warps or internally driven warps of the Galactic disk (Y. Revaz & D. Pfenniger 2004; J. A. Sellwood & V. P. Debattista 2022). However, none of these hypotheses have been quantitatively confirmed as the definitive origin of the MW's warp. Recent studies have suggested that the MW's S-shaped warp may be a long-standing structure, existing for over 5 Gyr (R. Roskar et al. 2010; M. López-Corredoira et al. 2014; H.-F. Wang et al. 2020; X. Li et al. 2023). Constructing a model of long-term evolution could provide insights into the origin of the Galactic warp, and such a period may coincide with the time that the Gaia–Sausage–Enceladus (GSE) merger completed.

Lines of evidence have shown that the GSE, which is considered the last major merger of our Galaxy, constitutes the

bulk of the inner halo (e.g., V. Belokurov et al. 2018; A. Helmi et al. 2018; R. P. Naidu et al. 2020). Consequently, the GSE has been studied qualitatively via analogs in cosmological simulations (L. A. Bignone et al. 2019; L. M. Elias et al. 2020), MW zoom-ins (A. Fattahi et al. 2019; R. J. J. Grand et al. 2020), and existing merger simulations (A. Helmi et al. 2018; H. H. Koppelman et al. 2020). These studies have achieved success in demonstrating how this merger produces the eccentric debris and reshapes the early MW disk.

Building upon these findings, R. P. Naidu et al. (2021) made a tailored model for this merger by producing a grid of 500 idealized galaxy merger simulations with the GADGET code (V. Springel 2005; V. Springel et al. 2021) to identify a fiducial model that best matches the H3 survey data (C. Conroy et al. 2019). They produced a configuration to replicate the H3 data and explain disparate phenomena across the Galaxy. Recent evidence has reported that the GSE is likely a significantly gas-rich merger. This was highlighted by I. Ciucă et al. (2024) in their analysis of the age–metallicity relationship from the APOGEE-2 DR17 survey (Abdurro'uf et al. 2022). They observed a notable hallmark of a gas-rich merger: an increase in [Mg/Fe] corresponding to a decrease in [Fe/H] around $\tau \approx 12$ Gyr, closely aligning with the GSE infall period. However, despite the success of the model proposed by R. P. Naidu et al. (2021), it did not fully account for the considerable amount of gas introduced into the Galaxy by the GSE merger, which significantly contributes to the growth of both the high-$\alpha$ and low-$\alpha$ disks.

Based on measurements of precession and analyses of differences in other warp characteristics using mono-age samples, some studies have proposed an externally excited formation mechanism for the disk warp (X. Cheng et al. 2020; E. Poggio et al. 2020). However, these studies generally prefer recent interactions with satellite galaxies that drive transient warps. A. Bosma (1991) found that at least half of spiral galaxies exhibit warps, suggesting a long-lived and universal warp mechanism. A hallmark feature of $\Lambda$CDM cosmology is a

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, People's Republic of China; ducuihua@ucas.ac.cn</footnote>

<footnote n="2">²GEPI, Observatoire de Paris, Universite PSL, CNRS, Place Jules Janssen, 92195 Meudon, France</footnote>
hierarchical assembly (e.g., S. D. M. White & C. S. Frenk 1991), which means the warp might be formed in a universal way by galaxy mergers.

In this study, we advance our understanding of the Galactic warp by creating a gas-rich GSE merger simulation to reconstruct the morphology of our Galaxy's disk. In Section 2, we provide a description of the simulation model and the initial conditions (ICs). In Section 3, we illustrate the long-lived and nonsteady warp model's evolution over time, comparing it with observational data at present. In Section 4, we calculate the warp precession rate and demonstrate its lopsided feature. In Section 5, we show that the disk is embedded in a live DM halo that is tilted and retrograde, possibly contributing to sustaining the long-lived warp. Finally, we discuss and summarize our results in Sections 6 and 7.

## 2. Model and Simulation

In R. P. Naidu et al. (2021), they made three different models to bracket the mass range of (2–7) × 10⁸ M☉ for the GSE progenitor. Each model considered three scale lengths to account for the significant scatter in size at fixed mass observed at z ~ 2. After simulation tests, they selected the fiducial model, determining the GSE progenitor mass to be 5 × 10⁸ M☉ and the scale length to be 1.5 × SMR (size–mass relation). The DM halo mass was derived from the z = 2 stellar mass-halo mass model in P. Behroozi et al. (2019). The MW was modeled as a combination of a thick disk and bulge with a total stellar mass of 2 × 10¹⁰ M☉, following the scale lengths in J. Bland-Hawthorn & O. Gerhard (2016), and the DM halo was set to half the z = 0 mass. More details of how these parameters are elaborated can be found in Section 4 in R. P. Naidu et al. (2021). In summary, they created faithful representations of both the GSE and the MW progenitors as follows: the GSE comprises a stellar disk and dark matter (DM) halo, while the MW comprises a bulge, stellar disk, and DM halo.

The numerical ICs are generated from the Disk Initial Conditions Environment (DICE) code (V. Perret et al. 2014; V. Perret 2016). The assumed density profiles for DM/stars/ gas are input into DICE as distribution functions, which then generate Lagrangian particles with the Metropolis–Hastings Markov Chain Monte Carlo algorithm (N. Metropolis et al. 1953). In our simulation, We add a gas disk to each progenitor. The MW progenitor comprises a bulge, a stellar disk, a gas disk, and a DM halo, whereas the GSE progenitor is modeled with a stellar disk, a gas disk, and a DM halo. Both disk models follow an exponential + sech-z profile. The bulge is modeled using an Einasto profile (J. Einasto 1965) to replace the Hernquist profile (L. Hernquist 1990) used in R. P. Naidu et al. (2021). The Einasto profile, which is similar to the Sérsic profile (J. L. Sérsic 1968) but used for 3D mass density (D. Coe 2010), is a reasonable choice since the Sérsic profile is widely used for fitting the surface density of galaxies. The DM halo is modeled using the Hernquist profile.

In Table 1, we present the ICs of our simulation, where the virial masses of the MW and GSE progenitors at z ~ 2 are set to 5 × 10¹¹ M☉ and 2 × 10¹¹ M☉, respectively. Considering the higher gas fractions in galaxies at z ~ 2, we set the GSE's gas fraction to 0.8 and the MW's to 0.5, slightly higher than the typical values observed at z ~ 1.5 (M. Rodrigues et al. 2012). We set the gas disk scale length for the MW and GSE progenitors to 3R<sub>s</sub> and 4R<sub>s</sub> (R<sub>s</sub> is the scale radius of the stellar disk), respectively. These choices follow model assumptions

<table><caption>**Table 1**
ICs at z ~ 2</caption><thead><tr><th>Progenitor Parameters</th><th>MW</th><th>GSE</th><th>Units</th></tr></thead><tbody><tr><td>DM mass</td><td>4.6</td><td>1.975</td><td>1 × 10<sup>11</sup> M<sub>☉</sub></td></tr><tr><td>DM scale radius</td><td>26.7</td><td>19.1</td><td>kpc</td></tr><tr><td>Stellar disk mass</td><td>60</td><td>5</td><td>1 × 10<sup>8</sup> M<sub>☉</sub></td></tr><tr><td>Stellar scale length</td><td>2</td><td>1.5</td><td>kpc</td></tr><tr><td>Stellar scale height</td><td>1</td><td>0.75</td><td>kpc</td></tr><tr><td>Gas disk mass</td><td>20</td><td>2</td><td>1 × 10<sup>9</sup> M<sub>☉</sub></td></tr><tr><td>Gas scale length</td><td>6</td><td>6</td><td>kpc</td></tr><tr><td>Gas scale height</td><td>3</td><td>3</td><td>kpc</td></tr><tr><td>Bulge mass</td><td>1.4</td><td>N/A</td><td>1 × 10<sup>10</sup> M<sub>☉</sub></td></tr><tr><td>Bulge scale</td><td>1.5</td><td>N/A</td><td>kpc</td></tr><tr><td colspan="4"><strong>Orbital Parameters</strong></td></tr><tr><td>Eccentricity (e)</td><td>0.75</td><td></td><td></td></tr><tr><td>Initial distance</td><td>141.2</td><td></td><td></td></tr><tr><td>Normal vector of orbit (θ, φ)</td><td>(0, 125)</td><td></td><td></td></tr><tr><td>MW and GSE spin (θ<sub>1</sub>, θ<sub>2</sub>, κ)</td><td>(125, 55, 180)</td><td></td><td></td></tr><tr><td>Relative velocity</td><td>73.02</td><td></td><td></td></tr></tbody></table>

made in T. J. Cox et al. (2008) and M. Rocha et al. (2008) that were used in simulation studies of galaxy mergers and metallicity gradients. The resolution mass for baryon particles is 1 × 10⁴ M☉ and for DM particles is 1 × 10⁵ M☉, and the initial relative velocity is calculated based on Keplerian orbit. The modeling approach utilized in this study is based on GIZMO (P. F. Hopkins 2015), in which we have implemented the module of star formation and feedback processes as in J. Wang et al. (2012), as star formation can convert gas particles into stellar particles and stars can turn back into gas via feedback. It is a variation of the smooth particle hydrodynamics code GADGET-3 and uses adaptive gravity softening (more details in P. F. Hopkins et al. 2018).

As these structural choices play a secondary role in the total mass and orbital parameters described in R. P. Naidu et al. (2021), they modeled the orbit with circularity: η = 0.5, inclination: θ = 15°. In V. Belokurov et al. (2023), they conducted a similar simulation broadly following R. P. Naidu et al. (2021), but changed the inclination to θ = 30°. Based on the fiducial model of R. P. Naidu et al. (2021), we aimed to reconstruct both the GSE and the warp structure in our simulation, leading to some differences in our orbital settings. Two galaxies with opposite disk spins were positioned in a radially biased and retrograde orbit. The initial distance between them was set equal to the sum of their virial radii. The eccentricity was set to 0.75, which is more radial than the moderate η = 0.5 orbit. The polar (θ) and azimuthal (φ) angles of the orbital plane's normal vector were set at (0°, 125°). This configuration indicates that the GSE is inclined at 35° relative to the MW disk plane in the Z-direction, slightly larger than found in V. Belokurov et al. (2023), and the spin angle is set to (θ₁, θ₂, κ) = (125°, 55°, 180°), θ₁ is the angle between the spin vector of the first galaxy and the orbital plane, θ₂ is the angle between the spin vector of the second galaxy and the orbital plane, and κ is the angle between the spin vector of the first galaxy and the second one (see Figure 2 in V. Perret et al. 2014 for more details), and the initial relative velocity is calculated with the orbital parameters. After the gas-rich merger happened, gas particles could inherit its angular momentum from the orbital momentum of the merger and then be redistributed into a thin disk (J. E. Barnes 2002). It will take
<figure><img src="image_3.png" /><figcaption>Figure 1. Distribution of simulated stars and gas in Galactocentric cylindrical coordinates at the present time (9.05 Gyr since the simulation started). Negative/positive R indicates the southern/northern warp. The top panel depicts the warp of the stellar disk, with a color bar indicating stellar ages. The bottom panel shows the gas warp in a density projection map with a contour profile. Here, the green lines represent a height of |Z| = 2.5 kpc, indicating the northern warp is stronger than the southern part. Comparative data are sourced from X. Chen et al. (2019), Lemasle (2022), and E. S. Levine et al. (2006).</figcaption></figure>

several Gyrs for the subsequent virialization phase to rebuild the disk (F. Hammer et al. 2009; P. F. Hopkins 2009), during which instabilities may develop some oscillations that we would like to search for. A more detailed summary of the constructed GSE properties is provided in the Appendix.

### 3. A Long-lived and Nonsteady Warp

The merger occurs rapidly, taking approximately 1 Gyr between the first and final pericenter passage, and is completed within the first 3 Gyr, which is shown in the orbital decay profile in panel of Figure 7 in the Appendix. The GSE merger is estimated to have occurred approximately about 8–11 Gyr ago (A. Helmi et al. 2018; V. Belokurov et al. 2018). In Figure 1, we provide an edge-on perspective of the Galactic disk, selecting particles based on their circularity: $\epsilon = L_Z/L_{Z,max(E)} > 0.65$. This selection criterion is applied at 9.05 Gyr after the simulation began, corresponding to what we expected as the present ($z \sim 0$). This time aligns with Auriga 18 (Au18), an MW-like disk galaxy from the Auriga project (R. J. J. Grand et al. 2017), which is frequently referenced in GSE research. This visualization distinctly reveals an S-shaped warp in both the stellar and gas disks, consistent with observational data. Stars are color-coded based on their ages; elder stars are represented with lighter green markers, while younger stars are depicted with more vibrant blue markers. The warp is the most pronounced for the youngest stars and spans multiple age groups (X. Chen et al. 2019; X. Cheng et al. 2020; Y. Huang et al. 2024). The geometric shape of the warp can be

approximated by a power-law warp mode:

$$Z_w(R \ge R_w) = a(R - R_w)^b \sin(\phi - \phi_w). \quad (1)$$

$$Z_w(R < R_w) = 0. \quad (2)$$

Here, R and $Z_w$ are Galactocentric cylindrical coordinates, a represents the amplitude of the warp, b is the power-law index, $\phi_w$ is the polar angle of the warp’s line of nodes (LONs), indicating the orientation of the warp, and $R_w$ is the onset radius of the warp. This model was employed by X. Chen et al. (2019) to accurately fit the Cepheid tracers, prompting us to select stars within the young stars’ age range ($age < 1$ Gyr) in the simulation for comparison with the Cepheid data. As X. Chen et al. (2019) note, the parameters of a, b, and $R_w$ show a clear correlation, whereas $\phi_w$ is a more independent parameter in the geometric warp model. Consequently, while different studies might yield varying a, b, and $R_w$ values, $\phi_w$ tends to be more consistent across measurements of the Galactic warp. In our analysis, we refrain from comparing $\phi_w$ with observational data, as its value can change based on the reference frame specified in the simulation, so we define the $\phi_w$ as the y-axis, making the orientation similar to that of the MW.

In the upper panel of Figure 1, we fitted the warp model in simulations and found that the linear model ($b=1$) closely matches the observational data; the simulation’s maximum warp amplitude is marked by a yellow dashed line, aligning closely with the Cepheid data (X. Chen et al. 2019; Lemasle 2022). The gas density projection map with a contour profile is displayed in the lower panel, in the central $\sim 1$ kpc where the formation of the Galactic bar creates a noticeable gap
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Evolution of stellar/gas disk warp amplitude as measured at R = 16 kpc, from 3 Gyr to 9.05 Gyr (as the merger occurred within the first 3 Gyr), and the stellar disk amplitude is also fitting with stars younger than 1 Gyr. Both disks exhibit similar evolution. Notably, since the disk had not fully shaped before 3.25 Gyr, there are two approximate amplitude peak values observed around 4.3 Gyr and 6.4 Gyr.</figcaption></figure>

in the inner region. Both sides of the gas warp align closely with the H I data from E. S. Levine et al. (2006) within a 15 kpc range. Beyond this region, observational data reveal that the northern section (with a positive peak) exhibits a higher amplitude, whereas the southern part (with a negative peak) displays a lower amplitude. Nonetheless, the solid line remains within the high gas density region. Despite these minor discrepancies, the overall results are deemed satisfactory. Previous research utilizing different tracer types has noted different observational systematics in the measured amplitudes (Z. He 2023).

Previous studies have reported the asymmetry of the warp in both gas and stellar components, with findings suggesting that the northern warp is more pronounced than its southern counterpart (E. S. Levine et al. 2006; C. Reylé et al. 2009; E. B. Amôres et al. 2017). Accordingly, we plot horizontal lines to represent the height of $|Z| = 2.5$ kpc in the bottom panel of Figure 1. This visualization shows the warp’s asymmetry in our simulation: the southern part falls slightly below this height, while the northern part is consistent with observational data. It should be noted that in our simulation, the disk warp evolves continuously, and such asymmetry could occasionally disappear or suddenly intensify.

Figure 2 illustrates the time evolution of the stellar/gas disk warp’s amplitude, calculated by fitting Equation (1) to the warp every 0.01 Gyr since $t=3$ Gyr after the merger has been completed and identifying the maximum amplitude at R = 16 kpc, a radius near the edge of the disk in our simulation, and the stellar disk is fit with stars younger than 1 Gyr. Both disks exhibit similar evolution, as shown in Figure 2, as gas particles can convert into young stars in the simulation, and such an evolution trend highlights the importance of the gas component in reconstructing the Galactic warp. The thin disk was gradually formed after the merger in our simulation, corresponding to a lookback time of about 6 Gyr, which is close to the observational time frame when the thin disk gradually formed about 7 Gyr ago (O. N. Snaith et al. 2014), the disk had not fully shaped before 3.25 Gyr, which led to chaotic conditions and potentially inaccurate measurements prior to this point. Afterward, the warp reached its initial peak at approximately 4.3 Gyr, decreased until around 5.2 Gyr, saw another peak at 6.4 Gyr, and then generally diminished to the present. Thus, we can identify roughly many evident extremum values in Figure 2, indicating that the warp is a long-lived and nonsteady structure.

## 4. Kinematic Warp Model

In Section 3, we discussed the evolution of the Galactic warp, highlighting how its rapid amplitude changes underscore the dynamical nature of warps. This dynamism is crucial for unlocking insights into the formation history of galaxies and the mass distribution of their halos. The changing geometry of the warp can be succinctly described by the variation in the direction of its LONs at a precession rate $\omega$. This is expressed as $\phi_w(t) = \phi_{w,0} + \omega t$, where $\phi_{w,0}$ represents the current position of the LONs. A time-dependent model of the warp is provided by E. Poggio et al. (2020), described as

$$ \bar{V}_Z(R, \phi, t=0) = \left( \frac{\bar{V}_\phi}{R} - \omega(R) \right) h_w(R) \cos(\phi - \phi_\omega) + \frac{\partial h_w}{\partial t} \sin(\phi - \phi_w), \quad (3) $$

where $h_w(R) = a(R - R_w)^b$, as described in Equation (1), $\bar{V}_Z$ represents the mean vertical velocity, and $\bar{V}_\phi$ is mean azimuthal velocity. To simplify the calculations, we neglect the time derivative term $\frac{\partial h_w}{\partial t}$ and assume $\omega$ does not vary with radius. Although the warp is evolving, our calculations indicate that this approximation does not significantly affect the precession results. With this assumption, the model could be described the same as Equation (7) in E. Poggio et al. (2020):

$$ \bar{V}_Z(R, \phi) = \left( \frac{\bar{V}_\phi}{R} - \omega \right) h_w(R) \cos(\phi - \phi_w). \quad (4) $$

With this kinematic model, after determining the geometric shape through fitting with Equation (1), we can plug in the parameters of $h_w$ and $\phi$ into Equation (4) and measure the precession rate $\omega$ using kinematic data from the simulation. Here, we select the young stars (age < 1 Gyr) within the radial range of R = 11~13 kpc and R = 13~15 kpc, plotting the evolution of $\omega$ over time in the top panel of Figure 3. The dashed line represents the nonprecession state, while positive and negative values of $\omega$ indicate prograde and retrograde warp precession relative to disk rotation, respectively. The trend of $\omega$ fluctuates between positive and negative values; for example, the R = 13~15 kpc precession trend ranges approximately from -13 to 18 km s⁻¹ kpc⁻¹. Although there are periods where the precession rate is notably high in either the prograde or retrograde direction, it typically remains within a range of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: the evolution of warp precession as measured with young stars within the radial range of *R* = 11 ~ 13 kpc and *R* = 13 ~ 15 kpc from 3 Gyr to 9.05 Gyr. Middle panel: the transparent band shows the distribution trend of the young stars' velocity in the Z-direction vs. azimuth within the region of *R* = 11 ~ 13 kpc (left) and *R* = 13 ~ 15 kpc (right); the purple line is the mean velocity, and the pink line indicates the fitting precession with kinematic warp model. Bottom panel: the vertical velocity differences of stars on both sides of *φ*<sub>w</sub> vs. azimuth angle separation from *φ*<sub>w</sub>. The azimuth angles are binned at 10°, and each radial annulus is plotted with different colors.</figcaption></figure>

−5 ~ 5 km s⁻¹ kpc⁻¹. A significant shift in the precession rate is observed around *t* = 5.5 Gyr in Figure 2, which can be attributed to the exceptionally low warp amplitude during this period. Comparing both the evolution trend of precession and amplitude in Figure 2, we could find that both consistently reach their extremum simultaneously, indicating a synchronous evolution. In the middle panel of Figure 3, the transparent band illustrates the velocity trend of simulated stars younger than 1 Gyr within the radial ranges of *R* = 11 ~ 13 kpc (left) and *R* = 13 ~ 15 kpc (right) at the present time, with the dashed line representing the mean velocity. This clearly demonstrates the warp's dynamical characteristics. In our simulation, it shows a high prograde precession in the range of *R* = 11 ~ 13 kpc. Incorporating all geometric model parameters into Equation (4), we obtain a precession rate of 13.820 and 4.795 km s⁻¹ kpc⁻¹ for two ranges; this value is close to some current precession observational data. It can be found that the fitting curve does not perfectly cover the mean velocity. This discrepancy arises because the vertical velocity is asymmetric about the longitude of peak vertical velocity, indicating that the warp is lopsided (E. B. Amôres et al. 2017; X. Cheng et al. 2020; X. Li et al.

2023). In the bottom panel of Figure 3 (refer to X. Cheng et al. 2020), we measure the velocity asymmetry by subtracting the median vertical velocity of stars on one side of $\phi_w$ from its complement on the other side at the same azimuthal separation in each radial annulus. A gradually increasing difference in vertical velocity is found when the azimuth angle is moving away from the velocity peak, consistent with the feature reported in X. Cheng et al. (2020). This trend may become more pronounced in the outer disk, as the disk is asymmetric. We also found that the lopsided feature is more evident in the $\phi < 0^\circ$ area.

In E. Poggio et al. (2020), precession was calculated using 12 million red giant stars from Gaia DR2 (Gaia Collaboration et al. 2018), yielding a precession rate of 10.86 ± 0.03<sub>stat</sub> ± 3.20<sub>syst</sub> km s<sup>-1</sup> kpc<sup>-1</sup> in the direction of Galactic rotation. Recently, a similar model was employed in X. Y. Zhou et al. (2024), using 134 Cepheids and combined with the radial velocity (RV) from Gaia DR3 (Gaia Collaboration et al. 2023), determining a lower warp precession rate of 4.9 ± 1.6 km s<sup>-1</sup> kpc<sup>-1</sup>. Our simulation is very close to this observation in different ranges. Another study by Ž.
<figure><img src="image_6.png" /><figcaption>Figure 4. Comparison between the obtained values of warp precession rate varying with Galactocentric radius R in our simulation (purple circle), together with other values from the literature.</figcaption></figure>

Chrobakova & M. Lopez-Corredoira (2021) also suggested a low precession rate and cannot exclude a nonprecessing warp of the MW. Consequently, determining precession accurately is hindered by the lack of comprehensive spatial and kinematic data, compounded by the need for reliable tracers. In Figure 4, we plot the variation of the precession rate with Galactocentric radius R. As shown, our simulation results align well with other studies in the literature (E. Poggio et al. 2020; M. Cabrera-Gadea et al. 2024; X. Y. Zhou et al. 2024) for R < ~15 kpc. Both W. Dehnen et al. (2023) and M. Cabrera-Gadea et al. (2024) reported a prograde precession that decreases with increasing Galactocentric radius, a trend that our simulation also captures. In our results, the precession rate continues to decrease, eventually turning negative (R ~ 15.5 kpc), indicating a transition to retrograde precession, as the outer disk is more susceptible to being disturbed and is also likely to be influenced by Sgr and LMC (I. A. Stelea et al. 2024), which might make the precession rate rise in observations beyond 15 kpc.

Theoretical and numerical studies of misaligned nonspherical halo potentials imply that the direction of precession can be influenced by the DM halo shape, as the prolate or prolate-like halo can rise prograde precession, while oblate-like halo should produce a retrograde warp precession (E. Poggio et al. 2020). However, warps in oblate halos are not appropriate to explain the observed warps or they could not sustain a prominent warp (M. Ideta et al. 2000; J. Dubinski & D. Chakrabarty 2009; M. Jeon et al. 2009). Typical precession periods from such models or those produced by a misaligned outer torus of later accreted material (I.-G. Jiang & J. Binney 1999; J. Shen & J. A. Sellwood 2006; M. Jeon et al. 2009) are approximately between 4 Gyr and 40 Gyr, corresponding to warp precession rates between 1.5 km s⁻¹ kpc⁻¹ and 0.1 km s⁻¹ kpc⁻¹. In Figure 3, the precession rate determined after the merger significantly surpasses all predicted values throughout the simulation period. Since the GSE merger led to an oblate halo that could not account for the observed characteristics of warp, a precession value greater than 10 km s⁻¹ kpc⁻¹ indicates a likely transient response observed in the outer disk. Such a response could be attributed to interactions with a satellite galaxy (M. D. Weinberg & L. Blitz 2006; C. F. P. Laporte et al. 2019), a scenario not included in our simulation. Given the warp's persistence over an extended period and its exhibition of both high prograde and retrograde precession post-merger, unraveling the underlying dynamical mechanisms poses a challenge.

## 5. Tilted and Retrograde DM halo

J. J. Han et al. (2022) discovered that the stellar halo, spanning a Galactocentric radius of 5–50 kpc, is tilted ~25° with respect to the disk plane. Following this discovery, they executed a numerical simulation utilizing a gravitational potential model of a galaxy. In their approach, 30% of the dark halo mass adopts a triaxial distribution, which is similarly tilted by 25° above the Galactic plane toward the Sun. Their results suggest that the warp and the flare of the Galactic disk might be related to the tilted halo (J. J. Han et al. 2023).

If so, our simulated disk might also be embedded in a tilted DM halo. We select the DM particles within 10 kpc < r < 50 kpc of the Galactocentric radius, and with the mass and position of each particle, we can calculate the moment of inertia tensor. Solving for the eigenvector-eigenvalue pairs of inertia tensor, we can find the three principal axes of rotation and their respective moments of inertia. The major axis has the minimum moment of inertia and the minor axis has the maximum moment of inertia. The length of the principal axes $r_i$, $i \in \{a, b, c\}$ has the relation with the moments of inertia $I_i$ as follows:

$$r_i^2 \propto \frac{-I_i + I_j + I_k}{2} \qquad (5)$$

Once we obtain the moments of inertia, we can calculate the major-to-intermediate and major-to-minor axes ratios as follows:

$$\frac{r_b}{r_a} : \frac{r_c}{r_a} = \sqrt{\frac{-I_b + I_a + I_c}{-I_a + I_b + I_c}} : \sqrt{\frac{-I_c + I_a + I_b}{-I_a + I_b + I_c}} \qquad (6)$$

Then, we can calculate the triaxiality with $T = (1 - p^2)/(1 - q^2)$, where p and q represent the ratios of major-to-intermediate and major-to-minor axes, respectively, $T > 0.6$ denotes a prolate halo, and $T < 0.3$ indicates an oblate halo. By measuring the orientation of the axes (with the disk plane's normal vector aligned parallel to the Z-axis), we are able to determine the tilt angle. Our results indicate that the major and intermediate axes are degenerate, a scenario also observed in S. Kazantzidis et al. (2004) and S. Shao et al. (2021). This degeneracy renders the measurement unstable if these axes are chosen. However, the minor axis is nondegenerate and can be reliably used to measure the tilt angle of the DM halo with respect to the disk plane.

In the middle and bottom panels of Figure 5, we illustrate the triaxiality with particles in 10 kpc < r < 50 kpc and the net angular momentum of the DM halo within 50 kpc. Following
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top panel: illustrating the evolution of the DM halo's tilt angle with respect to the disk. Middle panel: DM triaxiality variations with time, with values predominantly less than 0.3, indicating an oblate shape. Bottom panel: net angular momentum about the Z-axis $L_Z$ of the DM halo as a function of time. A negative value indicates a retrograde halo.</figcaption></figure>

the GSE merger, the simulation halo appears as an oblate profile, in alignment with G. Iorio & V. Belokurov (2019), and maintains a negative angular momentum, indicating a retrograde halo with respect to the disk. R. Joshi & L. M Widrow (2024) suggest that the vertical displacement of the disk in a retrograde halo will be stronger than in a prograde halo, leading to a more organized warp. A simple heuristic argument that dynamical friction tries to bring the angular momentum vectors of the halo and disk into alignment means the warp will be strongly damped in the corotating halo but enhanced in the counterrotating halo (R. W. Nelson & S. Tremaine 1995). Thus, the retrograde halo resulting from the GSE merger could sustain a long-lived warp.

The top panel of Figure 5 illustrates the time evolution of the DM halo tilt angle with respect to the disk, showing a continuously increasing trend, and there is no evident correlation with the amplitude evolution. In the simulation of J. J. Han et al. (2023), they fixed the halo's tilt angle at 25° and found that the warp needs 1.5 Gyr to reach a steady-state amplitude, and after that, there is no significant change. This suggests that a live halo could facilitate a fast and nonsteady evolution of the warp. The halo's tilt angle continues to increase following the merger, which, as suggested by J. J. Han et al. (2023), should result in an increasing warp amplitude. However, our simulation could produce a warp with a low amplitude at $t = 5.2$ Gyr and shows irregular evolution, as illustrated in Figure 2. Therefore, we posit that in galaxy models with a more realistic cosmological setting, the actual mechanisms driving the warp evolution are likely to be more intricate than suggested by idealized simulations based solely on gravitational potential.

Such changes in the tilted angle suggest that the DM halo might also exhibit rotational precession. We can approximate the DM halo as a rotating rigid body in the absence of external moments. In this scenario, we identify the body-fixed principal axes for the DM halo, resulting in a moment of inertia tensor in the simplified form of $[I] = \text{diag}(I_{xx}, I_{yy}, I_{zz})$. Consequently, the angular momentum vector can be expressed as $\mathbf{H}_G = \{0, 0, H_G\}$.

To describe the precession of the DM halo, we select the particles within 50 kpc and establish a body-fixed coordinate system by applying Euler angles with the rotation order $R_Z R_X R_Z$ to translate the X-Y-Z absolute reference frame system, which is initially defined in the simulation box. The free motions of the DM halo are then characterized by changes in these three Euler angles, which can be calculated as follows:

$$ \dot{\phi} = H_G \left( \frac{\cos^2 \psi}{I_{yy}} + \frac{\sin^2 \psi}{I_{xx}} \right), \quad (7) $$

$$ \dot{\theta} = H_G \left( \frac{1}{I_{xx}} - \frac{1}{I_{yy}} \right) \sin \theta \sin \psi \cos \psi, \quad (8) $$

$$ \dot{\psi} = H_G \left( \frac{1}{I_{zz}} - \frac{\cos^2 \psi}{I_{yy}} - \frac{\sin^2 \psi}{I_{xx}} \right) \cos \theta. \quad (9) $$

The $\psi, \phi, \theta$ are described as “roll,” “yaw,” and “pitch” angle, and $\dot{\psi}, \dot{\phi}, \dot{\theta}$ are spin, precession, and nutation rate, respectively. Under the condition that the DM halo exhibits an oblate profile, we align the minor axis with the z-axis in the body-fixed coordinate system, the evolution of three velocities are shown in Figure 6. The precession rate of the DM halo is along the original Z-axis. Given that the angular momentum of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Three velocities vary with time in different panels. The precession rate is retrograde relative to the disk. The spin velocity is along the body-fixed z-axis. The nutation velocity is along the x-axis after the first rotating with the *φ* angle.</figcaption></figure>

DM halo is opposite to that of the disk, the DM halo undergoes retrograde precession relative to the disk. This interaction may cause the disk warp to also exhibit retrograde precession, as shown in Figure 3, where the precession remains predominantly retrograde throughout most of the period.

## 6. Discussion

It is crucial to note that in our simulation, one of the important parameters influencing warp creation is the azimuthal angle of the normal vector of the orbital plane. Specifically, we chose a slightly larger angle of 125°, which resulted in a particularly small part of the GSE debris following a high-energy orbit. These stars do not align well with current observations, suggesting that the simulation configuration may not be the most accurate model for reconstructing the GSE. Additionally, our work does not consider other potential mechanisms for the creation of warp, such as interactions with satellite galaxies, or other minor/mini infall events. Therefore, this suggests that a single major merger model cannot simultaneously fit both the GSE and the disk warp very accurately.

In Section 4, we calculated the precession rate versus time, revealing the precession can alternate between prograde and retrograde. In the present, it aligns well with observational data for $R < \sim 15$ kpc. However, beyond this range our simulation shows a transition to retrograde precession, potentially inconsistent with current observational data in this range. However, precise measurements necessitate extensive kinematic data. Additionally, gravitational interactions with a satellite galaxy, not accounted for in our simulation, could induce a prograde precession. Given the extensive evolutionary history of the MW, it is plausible that other mechanisms

affecting the precession rate have occurred. Therefore, incorporating additional models partway through this period may provide a more comprehensive understanding of these dynamics.

In Section 5, we find our disk is embedded in a time-evolving tilted and retrograde DM halo. In J. J. Han et al. (2023), their research focused exclusively on a fixed tilt angle of the DM halo, neglecting the long-term dynamical mechanisms that emerge from mergers. As they suggest that the warp requires 1.5 Gyr to reach a steady state, our simulations indicate that such an extended period could result in significant changes within the DM halo. Therefore, we posit that variations of the tilt angle relative to the disk might influence the formation of warp. Notably, the time evolution of the DM halo tilt angle does not exhibit a clear correlation with the observed changes in warp amplitude. This implies that in a comprehensive galaxy model, the impact of a tilted DM halo on warp formation may not be as straightforward as predicted by idealized gravitational potential models. Consequently, the Galactic warp model appears to necessitate a series of controlled experiments for further investigation, which we are preparing to conduct in our next study.

## 7. Conclusion

We conducted a hydrodynamical simulation of a gas-rich GSE merger, positioned on reconstructing the GSE debris, successfully producing the Galactic warp amplitude and precession that aligns well with observational data, discovering that the warp is long-lived, nonsteady, asymmetric, and lopsided, exhibiting both a prograde and retrograde precession rate after the completion of the merger. We found significant changes in the DM halo as a consequence of the major merger,
which created a live, oblate, tilted, and retrograde DM halo. Specifically, the tilted and retrograde characteristics of the DM halo appear to sustain the warp.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NSFC Nos. 11973042 and 11973052). We thank Xiaodian Chen for sharing the data. We are grateful to Phil Hopkins and Jianling Wang, who kindly shared with us access to the Gizmo code. We are grateful for the support of the International Research Program Tianguan, which is an agreement between the CNRS in France, NAOC, IHEP, and Yunnan University in China.

Appendix

Characteristics of the Simulated GSE

Here, we list a summary of our GSE model in Figure 7. While most of the galaxy model parameters are adapted from R. P. Naidu et al. (2021), our model features a significantly different orbit configuration. Nevertheless, our model successfully reconstructs the global characteristics of the GSE debris. The top-left panel illustrates the orbital decay profile. Owing to the large eccentricity of the orbit setting, the merger occurs rapidly, taking approximately only 1 Gyr between the first and final pericenter passages. The top-middle panel shows the all-sky density profile of the GSE, and the power-law coefficients are given by the $\rho \propto r_{gal}^{\alpha}$ profile. As R. P. Naidu et al. (2021) proposed, we also adopted a double-break density profile, which is used widely for the inner halo ($r_{gal} < 30$ kpc), and

noticed that there are two breaks in $r_{gal} = 15$ kpc and $r_{gal} = 30$ kpc, respectively. The break at 15 kpc is a critical feature of the GSE merger (R. P. Naidu et al. 2021) and it is also roughly the radius where the disk warp reaches maximum amplitude, which may hint at a correlation between the GSE and warp. Since the GSE is by far the most dominant component of the inner halo, we expect the overall halo density profile at this range could be fitted in our simulation. We find the slope of the GSE density profile between two breaks is a good match to that found for the inner halo (A. J. Deason et al. 2014), and the third coefficient after the second break fits well with that of D. Ye et al. (2023), as they report the second break at 27.18 kpc with $\alpha = 2.86$. Meanwhile, the top-right panel displays the anisotropy parameter $\beta$ of the GSE debris with Galactocentric distance in a spherical coordinate system, which is defined as

$$ \beta = 1 - \frac{\sigma_{\theta}^{2} + \sigma_{\phi}^{2}}{2\sigma_{r}^{2}}, \qquad (A1) $$

where $\beta$ quantifies the degree of velocity anisotropy of a system of stellar orbits, where $\sigma_i$ is the velocity dispersion in spherical coordinates. The observational data is from W. B. Wu et al. (2022), derived from the LAMOST (G. Zhao et al. 2006, 2012; X.-Q. Cui et al. 2012; A. L. Luo et al. 2012) K-giant Sgr-removed stellar halo. There exists a disparity between the parameter values in the data and the simulation results. However, both exhibit a consistent trend characterized by an increasing anisotropy parameter within 17 kpc and a

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top-left panel: the orbital decay profile, the merger completed within the first 3 Gyr. Top-middle panel: the all-sky density profile of the GSE (blue points/lines). There are two breaks in our simulation and the red dashed line is the observational stellar halo profile (A. J. Deason et al. 2014). Top-right panel: the anisotropy parameter $\beta$ of the GSE debris with the Galactocentric distance in a spherical coordinate system, compared with observational data. Bottom panels: the most representative characteristic of the GSE potted with a contour profile. The “sausage” structure is in the left panel and the distribution of the particle binding energy vs. angular momentum is in the right panel. Here, the negative value of $V_\phi$ and $L_z$ means that the stars on a retrograde orbit. Except for the top-left panel, all the other panels are plotted at 9.05 Gyr.</figcaption></figure>
relatively stable to decreasing trend between approximately 17 and 28 kpc, while increasing after that. The bottom panels show the most prominent features of the GSE in both the $V_φ - V_r$ panel and the $E - L_z$ panel. We can find the “sausage” structure in the velocity panel, which was first discovered by V. Belokurov et al. (2018), and is where the “Gaia-Sausage” came from. The contour profile tends to have negative values (indicating retrograde orbits in our setting) on both sides, implying a slight overall retrograde motion. The same feature is also reflected in the $E - L_z$ panel, where the majority of stars are situated in the central region, and the profile similarly trends toward negative values.

## ORCID iDs

Mingji Deng https://orcid.org/0000-0002-0592-7660

Cuihua Du https://orcid.org/0000-0002-3954-617X

Yanbin Yang https://orcid.org/0000-0001-7949-3407

Jiwei Liao https://orcid.org/0000-0001-6762-5599

## References

Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35

Amôres, E. B., Robin, A. C., & Reylé, C. 2017, A&A, 602, A67

Bailin, J. 2003, ApJ, 583, L79

Barnes, J. E. 2002, MNRAS, 333, 481

Battaner, E., & Jiménez-Vicente, J. 1998, A&A, 332, 809

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Belokurov, V., Vasiliev, E., Deason, A. J., et al. 2023, MNRAS, 518, 6200

Bignone, L. A., Helmi, A., & Tissera, P. B. 2019, ApJ, 883, L5

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bosma, A. 1991, Warped and Flaring HI Disks, 181 (Cambridge: Cambridge Univ. Press)

Chrobáková, Ž., & López-Corredoira, M. 2021, ApJ, 912, 130

Cabrera-Gadea, M., Mateu, C., Ramos, P., et al. 2024, MNRAS, 528, 4409

Chen, X., Wang, S., Deng, L., et al. 2019, NatAs, 3, 320

Cheng, X., Anguiano, B., Majewski, S. R., et al. 2020, ApJ, 905, 49

Ciucă, I., Kawata, D., Ting, Y.-S., et al. 2024, MNRAS, 528, L122

Coe, D. 2010, arXiv:1005.0411

Conroy, C., Bonaca, A., Cargile, P., et al. 2019, ApJ, 883, 107

Cox, T. J., Jonsson, P., Somerville, R. S., Primack, J. R., & Dekel, A. 2008, MNRAS, 384, 386

Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, RAA, 12, 1197

Deason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30

Dehnen, W., Semczuk, M., & Schönrich, R. 2023, MNRAS, 523, 1556

Dubinski, J., & Chakrabarty, D. 2009, AJ, 703, 2068

Einasto, J. 1965, TrAlm, 51, 87

Elias, L. M., Sales, L. V., Helmi, A., & Hernquist, L. 2020, MNRAS, 495, 29

Fattahi, A., Belokurov, V., Deason, A. J., et al. 2019, MNRAS, 484, 4471

Freudenreich, H. T., Berriman, G. B., Dwek, E., et al. 1994, ApJL, 429, L69

Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 616, A11

Gaia Collaboration, Drimmel, R., Romero-Gómez, M., et al. 2023, A&A, 674, A37

Grand, R. J. J., Gómez, F., Marinacci, F., et al. 2017, MNRAS, 467, 179

Grand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603

Hammer, F., Flores, H., Puech, M., et al. 2009, A&A, 507, 1313

Han, J. J., Conroy, C., & Hernquist, L. 2023, NatAs, 7, 1481H

Han, J. J., Conroy, C., Johnson, B. D., et al. 2022, AJ, 164, 249

He, Z. 2023, ApJ, 954L, 9H

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Hernquist, L. 1990, ApJ, 356, 359

Hopkins, P. F. 2009, in ASP Conf. Ser. Vol. 419, Galaxy Evolution: Emerging Insights and Future Challenges. Astron. Soc. Pac., San Francisco, ed. S. Jogee et al. (San Francisco, CA: ASP), 228

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., Wetzel, A., & Kereš, D. 2018, MNRAS, 480, 800

Huang, Y., et al. 2024, NatAs, in press

Ideta, M., Hozumi, S., Tsuchiya, T., & Takizawa, M. 2000, MNRAS, 311, 733

Iorio, G., & Belokurov, V. 2019, MNRAS, 482, 3868

Jeon, M., Kim, S. S., & Ann, H. B. 2009, ApJ, 696, 1899

Jiang, I.-G., & Binney, J. 1999, MNRAS, 303, L7

Joshi, R., & Widrow, L. M. 2024, MNRAS, 527, 7781J

Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJL, 611, L73

Kerr, F. J. 1957, AJ, 62, 93

Koppelman, H. H., Bos, R. O. Y., & Helmi, A. 2020, A&A, 642, L18

Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134

Lemasle, B., Lala, H. N., Kovtyukh, V., et al. 2022, A&A, 668, A40

Levine, E. S., Blitz, L., & Heiles, C. 2006, Sci, 312, 1773

Li, X., Wang, H.-F., Luo, Y.-P., et al. 2023, ApJ, 943, 88

López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101

López-Corredoira, M., Betancort-Rijo, J., & Beckman, J. E. 2002, A&A, 386, 169

Luo, A. L., Zhang, H.-T., Zhao, Y.-H., et al. 2012, RAA, 12, 1243

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, JChPh, 21, 1087

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2021, ApJ, 923, 92

Nelson, R. W., & Tremaine, S. 1995, MNRAS, 275, 897

Ostriker, E. C., & Binney, J. J. 1989, MNRAS, 237, 785

Perret, V., 2016 DICE: Disk Initial Conditions Environment, Astrophysics Source Code Library, ascl:1607.002

Perret, V., Renaud, F., Epinat, B., et al. 2014, A&A, 562, A1

Poggio, E., Drimmel, R., Andrae, R., et al. 2020, NatAs, 4, 590

Quinn, T., & Binney, J. 1992, MNRAS, 255, 729

Reshetnikov, V., & Combes, F. 1998, A&A, 337, 9

Revaz, Y., & Pfenniger, D. 2004, A&A, 425, 67

Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819

Rocha, M., Jonsson, P., Primack, J. R., & Cox, T. J. 2008, MNRAS, 383, 1281

Rodrigues, M., Puech, M., Hammer, F., Rothberg, B., & Flores, H. 2012, MNRAS, 421, 2888

Roskar, R., Debattista, V. P., Brooks, A. M., et al. 2010, MNRAS, 408, 783

Sánchez-Saavedra, M. L., Battaner, E., & Florido, E. 1990, MNRAS, 246, 458

Sánchez-Saavedra, M. L., Battaner, E., Guijarro, A., López-Corredoira, M., & Castro-Rodríguez, N. 2003, A&A, 399, 457

Sellwood, J. A., & Debattista, V. P. 2022, MNRAS, 510, 1375

Sérsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico)

Shao, S., Cautun, M., Deason, A., & Frenk, C. S. 2021, MNRAS, 504, 6033

Shen, J., & Sellwood, J. A. 2006, MNRAS, 370, 2

Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871

Stelea, I. A., Hunt, J. A. S., & Johnston, K. V. 2024, arXiv:2407.12095

Wang, H.-F., Lopez-Corredoira, M., & Huang, Y. 2020, ApJ, 897, 119

Wang, J., Hammer, F., Athanassoula, E., et al. 2012, A&A, 538, A121

Weinberg, M. D., & Blitz, L. 2006, ApJL, 641, L33

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

Wu, W. B., Zhao, G., Xue, X. X., et al. 2022, ApJ, 924, 23W

Ye, D., Du, C., Shi, J., & Ma, J. 2023, MNRAS, 525, 2472

Zhao, G., Chen, Y.-Q., Shi, J.-R., et al. 2006, ChJAA, 6, 265

Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, RAA, 12, 723

Zhou, X. Y., Chen, X. D., Deng, L. C., et al. 2024, ApJ, 965, 132