Previous studies have revealed that the Galactic warp is a longlived nonsteady and asymmetric structure. There is a need for a model that accounts for the warps longterm evolution. Given that this structure has persisted for over 5 Gyr its timeline may coincide with the completion of the GaiaSausageEnceladus GSE merger. Recent studies indicate that the GSE the significant merger of our Galaxy was likely a gasrich 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 gasrich 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 warps longterm evolution may offer more clues into the formation history of the Milky Way.
['Galaxy structure', 'Milky Way disk', 'Milky Way dark matter halo', 'Galaxy mergers', 'Hydrodynamical simulations', '622', '1050', '1049', '608', '767', 'Astrophysics - Astrophysics of Galaxies']
A Potential Dynamical Origin of the Galactic Disk Warp The GaiaSausageEnceladus Major Merger
{'A Potential Dynamical Origin of The Galactic Disk Warp: The Gaia-Sausage-Enceladus Major Merger': 'Mingji Deng, 1 Cuihua Du, 1 Yanbin Yang, 2 Jiwei Liao, 1 and Dashuang Ye 1  \n1 School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, P.R. China 2 GEPI, Observatoire de Paris, Universite PSL, CNRS, Place Jules Janssen, 92195 Meudon, France', '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 Gyrs, its timeline may coincide with the completion of Gaia-SausageEnceladus (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 that aligns closely with observational data from both stellar and gas disks. This simulation demonstrates the possibility that the 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.  \nKeywords: 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 (S'anchez-Saavedra et al. 1990, 2003; Reshetnikov & Combes 1998). The Milky Way, a typical example of disk and spiral galaxies, also exhibits a clear disk warp, confirmed by various works (Kerr 1957; Levine et al. 2006; Freudenreich et al. 1994; Chen et al. 2019). Many mechanisms of the Galactic warp have been proposed, including: inflow of intergalactic matter into the halo (Ostriker & Binney 1989; Quinn & Binney 1992; Jiang & Binney 1999); inflow directly onto the Galactic disk (L'opez-Corredoira et al. 2002a); magnetic fields that exist between galaxies (Battaner & Jim'enez-Vicente 1998); interaction between satellite galaxies like Sagittarius (Bailin 2003) or Magellanic Clouds (Weinberg & Blitz 2006) with the disk; and the bending instability and self-excited warps or internally driven warps of the Galactic disk (Revaz & Pfenniger 2004; Sellwood & Debattista 2022). However, none of these hypotheses have been quantitatively confirmed as the definitive origin of the Milky Way's warp. Recent studies have suggested that the MW's S-shaped warp may be a long-standing structure, existing for over 5 Gyrs (Roskar et al. 2010; L'opez-Corredoira et al. 2014; 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 Gaia-Sausage-Enceladus (GSE) merger completed.  \nLines of evidence have shown that the GSE, which is considered as the last major merger of our Galaxy, constitutes the bulk of the inner halo (e.g. Belokurov et al. 2018; Helmi et al. 2018; Naidu et al. 2020). Consequently, the GSE has been studied qualitatively via analogs in cosmological simulations (Bignone et al. 2019; Elias et al. 2020), MW zoom-ins (Fattahi et al. 2019; Grand et al. 2020), and existing merger simulations (Helmi et al. 2018; Koppelman et  \nTable 1. Initial Conditions at z ∼ 2  \nal. 2020). These studies have achieved success in demonstrating how this merger produces the eccentric debris and reshapes the early Milky Way disk.  \nBuilding upon these findings, 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 (Springel et al. 2005, 2021) to identify a fiducial model that best matches the H3 survey data (Conroy et al. 2019a). 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 Ciuc˘a 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 with a decrease in [Fe/H] around τ ≈ 12 Gyrs, closely aligning with GSE infall period. However, despite the success of the model proposed by 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α and lowα disks.  \nBased 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 (Poggio et al. 2020; Cheng et al. 2020). However, these studies generally prefer a recent interactions with satellite galaxies that drive transient warps. Bosma (1991) found that at least half of spiral galaxies exhibit warps, suggesting a long-lived and universal warp mechanism. A hallmark feature of ΛCDM cosmology is hierarchical assembly (e.g. White & Frenk 1991), which means the warp might be formed in a universal way by galaxies merger.  \nIn 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. 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 sustain the long-lived warp. Finally, we discuss and summarize our results in Section 6 and Section 7.", '2. MODEL AND SIMULATION': "In Naidu et al. (2021), they made three different models to bracket mass range of (2 -7) × 10 8 M ⊙ for the GSE progenitor. Each model considered three scale lengths to account for the significant scatter in size at fixed mass  \nobserved at z ∼ 2. After simulation tests, they selected the fiducial model, determining the GSE progenitor mass to be 5 × 10 8 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 Behroozi et al. (2019). The Milky Way was modeled as a combination of a thick disk and bulge with a total stellar mass of 2 × 10 10 M ⊙ , following the scale lengths in Bland-Hawthorn & Gerhard (2016), and the DM halo was set to half the z = 0 mass. More details of how these parameters are elaborated could be found in Section 4 in Naidu et al. (2021). In summary, they created faithful representations of both the GSE and the Milky Way progenitors as follows: the GSE comprises a stellar disk and dark matter halo, while the Milky Way comprises a bulge, stellar disk, and dark matter halo.  \nThe numerical initial conditions (ICs) are generated from the disk initial conditions environment (DICE) code (Perret et al. 2014; 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 Monte Carlo Markov Chain algorithm (Metropolis et al. 1953). In our simulation, We add a gas disk to each progenitors. The Milky Way progenitor comprises a bulge, a stellar disk, a gas disk, and a DM halo; where as 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 (Einasto 1965) to replace the Hernquist profile (Hernquist 1990) used in Naidu et al. (2021). The Einasto profile, which is similar to the S'ersic profile (S'ersic 1968) but used for 3-D mass density (Coe 2010), is a reasonable choice since the S'ersic profile is widely used for fitting the surface density of galaxies. The DM halo is modeled using the Hernquist profile.  \nIn Table 1, we present the initial conditions of our simulation, where the virial masses of the Milky Way and GSE progenitors at z ∼ 2 are set to 5 × 10 11 M ⊙ and 2 × 10 11 M ⊙ , respectively. Considering the higher gas fractions in galaxies at z ∼ 2, we set the GSE's gas fraction to 0.8 and the Milky Way's to 0.5, slightly higher than the typical values observed at z ∼ 1 . 5 (Rodrigues et al. 2012). We set the gas disk scale length for the Milky Way and GSE progenitors to 3 R s and 4 R s ( R s is the scale radius of the stellar disk), respectively. These choices follow model assumptions made in Cox et al. (2008) and 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 4 M ⊙ and for DM particles is 1 × 10 5 M ⊙ , and the initial relative velocity is calculated base on Keplerian orbit. The modeling approach utilized in this study is based on GIZMO (Hopkins 2015), in which we have implemented the module of star formation and feedback processes as in 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 Hopkins et al. (2018)).  \nAs these structural choices play a secondary role to the total mass and orbital parameters described in Naidu et al. (2021), they modeled the orbit with circularity: η = 0 . 5, inclination: θ = 15 · . In Belokurov et al. (2023), they conducted a similar simulation broadly following Naidu et al. (2021), but changed the inclination to θ = 30 · . Based on the fiducial model of Naidu et al. (2021), we aimed to reconstruct both the GSE and the warp structure in our the simulation, leading to some differences in our orbital settings. Two galaxies with opposite disk spins were positioned in a 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 GSE is inclined at 35 · relative to the MW disk plane in Z direction, slightly larger than Belokurov et al. (2023), and the spin angle is set to ( θ 1 , θ 2 , κ ) = (125 · , 55 · , 180 · ), θ 1 is the angle between the spin vector of the first galaxy and the orbital plane, θ 2 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 Fig.2 in Perret et al. (2014) for more details). 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 (Barnes 2002). It will take several Gyrs for subsequent virialization phase to rebuild the disk (Hopkins et al. 2009; Hammer et al. 2009), during which instabilities may develop some oscillations that we would like to search for. A more detailed summary of the constructed the GSE properties is provided in Appendix A.", '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 Gyrs, which is shown in the orbital decay profile panel of Figure 7. As the GSE merger is estimated to have occurred approximately about 8-11 Gyrs ago (Belokurov et al. 2018; Helmi et al. 2018).  \nFigure 1. Distribution of simulated stars and gas in Galactocentric cylindrical coordinates at the present time (9.05 Gyrs 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 density projection map with 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 Chen et al. (2019), Lemasle et al. (2022) and Levine et al. (2006).  \n<!-- image -->  \nFigure 2. Evolution of stellar/gas disk warp amplitude as measured at R = 16kpc, from 3 Gyr to 9.05 Gyr (as the merger occurred within the first 3 Gyrs), 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.  \n<!-- image -->  \nIn Figure 1, we provide an edge-on perspective of the galactic disk, selecting particles based on their circularity: ϵ = L Z /L Z,max ( E ) > 0 . 65. This selection criterion is applied at 9.05 Gyrs after simulation began, corresponding to what we expected as the present ( z ∼ 0). This time aligns with Auriga 18 (Au18), a Milky Way-like disk galaxy from the Auriga project (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  \nbased 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 (Chen et al. 2019; Cheng et al. 2020; Huang et al. 2024). The geometric shape of the warp can be approximated by a power-law warp mode:  \nZ w ( R ≥ R w ) = a ( R -R w ) b sin ( ϕ -ϕ w ) . (1)  \nZ w ( R < R w ) = 0 . (2)  \nHere, R and Z w are Galactocentric cylindrical coordinates, a represents the amplitude of the warp, b is the power-law index, ϕ w is the polar angel of warp's lines of node (LON hereafter) indicates the orientation of the warp and R w is the onset radius of the warp. This model was employed by 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 Chen et al. (2019) notes, the parameters of a , b and R w show a clear correlation, whereas ϕ w is a more independent parameter in the geometric warp model. Consequently, while different studies might yield varying a , b and R w values, ϕ w tends to be more consistent across measurements of the Galactic warp. In our analysis, we refrain from comparing ϕ w with observational data, as its value can change based on the reference frame specified in the simulation, so we define the ϕ w as the y-axis making the orientation similar to that of the MW.  \nIn 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 (Chen et al. 2019; Lemasle et al. 2022). The gas density projection map with contour profile is displayed in the lower panel, in central ∼ 1 kpc where the formation of the Galactic bar creates a noticeable gap in the inner region. Both sides of the gas warp align closely with the HI data from Levine et al. (2006) within a 15 kpc range. Beyond this region, observational data reveal that the northern section (with positive peak) exhibits a higher amplitude, whereas the southern part (with negative peak) displays a lower amplitude. Nonetheless, the solid line remains within 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 (He 2023).  \nPrevious 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 (Levine et al. 2006; Reyl'e et al. 2009; Amˆores 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 noticed that in our simulation, the disk warp evolves continuously, such asymmetry could occasionally disappear or suddenly intensify.  \nFigure 2 illustrates the time evolution of the stellar/gas disks warp's amplitude, calculated by fitting Eq. 1 to the warp every 0.01 Gyr since t = 3Gyr that after the merger has been completed and identifying the maximum amplitude at R = 16kpc, a radius near the edge of the disk in our simulation, and the stellar disk is fitting with stars younger than 1 Gyr. Both disks exhibit similar evolution in Figure 2, as gas particles can convert into young stars in the simulation, such evolution trend highlights the importance of gas component in reconstructing the Galactic warp. The thin disk was gradually formed after the merger in our simulation, corresponding to a look-back time of about 6 Gyr, which is close to the observational time frame when the thin disk gradually formed about 7 Gyr ago (Snaith et al. 2014), disk had not fully shaped before 3.25 Gyr, which leading to chaotic conditions and potentially inaccurate measurements prior to this point. Afterwards, 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 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 discuss 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 ω . This is expressed as ϕ w ( t ) = ϕ w, 0 + ωt , where ϕ w, 0 represents the current position of the LONs. A time dependent model of warp is provided by Poggio et al. (2020), described as:  \nV Z ( R,ϕ,t = 0) = ( V ϕ R -ω ( R )) h w ( R ) cos ( ϕ -ϕ ω ) + ∂h w ∂t sin ( ϕ -ϕ w ) . (3)  \nFigure 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 Z direction versus 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 ϕ w vs. azimuth angle separation from ϕ w . The azimuth angle are binned in 10 · , and each radial annulus are plotted with different colors.  \n<!-- image -->  \nWhere h w ( R ) = a ( R -R w ) b , as described in Eq. 1, V Z represents mean vertical velocity and V ϕ is mean azimuthal velocity. To simplify the calculations, we neglect the time derivative term ∂h w ∂t and assume ω 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 Poggio et al. (2020):  \nV Z ( R,ϕ ) = ( V ϕ R -ω ) h w ( R ) cos ( ϕ -ϕ w ) . (4)  \nWith this kinematic model, after determining the geometric shape through fitting with Eq. 1, we can plug in the parameters of h w and ϕ into Eq. 4 and measure the procession rate ω 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 ω over time in the top panel of Figure 3. The dashed line represents the non-precession state, while positive and negative values of ω indicate prograde and retrograde warp precession relative to disk rotation, respectively. The trend of ω fluctuates between positive and negative values, for example the R = 13 ∼ 15 kpc precession trend ranges approximately from -13 to 18 km s -1 kpc -1 . 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 -5 ∼ 5kms -1 kpc -1 . A significant shift in the precession rate is observed around t = 5 . 5 Gyr in Figure 2, which can be attributed to the exceptionally  \nFigure 4. Comparison between the obtained values of warp precession rate varies with Galactocentric radius R in our simulation (purple circle), together with other values from the literature.  \n<!-- image -->  \nlow warp amplitude during this period. Comparing both the evolution trend of precession and amplitude in Figure 2, we could find that both two 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 Eq. 4, we obtain a precession rate of 13.820 and 4.795 km s -1 kpc -1 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 (Amˆores et al. 2017; Cheng et al. 2020; Li et al. 2023). In the bottom panel of Figure 3, refer to Cheng et al. (2020) we measure the velocity asymmetry by subtracting the median vertical velocity of stars on one side of ϕ 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 moving away from the velocity peak, consistent with the feature reported in Cheng et al. (2020). This trend may become more pronounced in the outer disk, as the disk is asymmetric we can also found the lopsided feature is more evident in the ϕ < 0 · area.  \nIn 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 stat ± 3 . 20 syst kms -1 kpc -1 in the direction of Galactic rotation. Recently, the similar model was employed in Zhou et al. (2024), using 134 Cepheids and combined with the lineof-sight velocity (RV) from Gaia DR3 (Gaia Collaboration et al. 2023a), determining a lower warp precession rate of 4 . 9 ± 1 . 6 kms -1 kpc -1 , our simulation is very close to these observation in different ranges. Another study by Chrob'akov'a et al. (2021) also suggested a low precession rate and cannot exclude a non-precessing warp of 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 (Poggio et al. 2020; Zhou et al. 2024; Cabrera-Gadea et al. 2024) for R < ∼ 15 kpc. Both Dehnen et al. (2023) and 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 be disturbed, which is also likely to be influenced by Sgr. and LMC (Stelea et al. 2024), which might make precession rate rise in observation beyond 15 kpc.  \nTheoretical and numerical studies of misaligned non-spherical 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 oblatelike halo should produce a retrograde warp precession 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 (Dubinski, J. & Chakrabarty 2009; Jeon et al. 2009; Ideta et al. 2000). Typical precession periods from such models or those produced by a misaligned outer torus of later accreted material (Shen & Sellwood 2006; Jiang & Binney 1999; Jeon et al. 2009), are approximately between 4 Gyr and 40 Gyr, corresponding to warp precession rates between 1 . 5 kms -1 kpc -1 and 0 . 1 kms -1 kpc -1 . In Figure 3, the precession rate determined after the merger significantly surpasses all predicted values throughout the  \nFigure 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. Negative value indicates a retrograde halo.  \n<!-- image -->  \nsimulation period. Since GSE merger led to a oblate halo which could not account for the observed characteristics of warp, and a percession value greater than 10 km s -1 kpc -1 indicates a likely transient response observed in the outer disk. Such a response could be attributed to interactions with a satellite galaxy (Weinberg & Blitz 2006; 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': "Han et al. (2022) discovered that the stellar halo, spanning a Galactocentric radius of 5-50 kpc, is tilted approximately ∼ 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% the dark halo mass adopts a triaxial distribution, which is similarly tilted by 25 · above the Galactic plane towards the Sun. Their results suggest that the warp and the flare of the Galactic disk might be related to the tilted halo (Han et al. 2023).  \nIf 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 Galactocentric radius, with the mass and position of each particles 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 minor axis has the maximum moment of inertia. The length of the principal axes r i , i ∈ { a, b, c } has the relation with the moments of inertia I i as following (Peraire & Widnall 2009):  \nr 2 i ∝ -I i + I j + I k 2 . (5)  \nOnce we obtain the moments of inertia, we can calculate the major-to-intermediate and major-to-minor axes ratios as follows:  \nr b r a : r c r a = √ -I b + I a + I c -I a + I b + I c : √ -I c + I a + I b -I a + I b + I c . (6)  \nFigure 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, X axis after the first rotating with the ϕ angle.  \n<!-- image -->  \nThen we can calculate the triaxiality with T = (1 -p 2 ) / (1 -q 2 ), where p and q represent the ratios of major-tointermediate 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 Kazantzidis et al. (2004) and Shao et al. (2021). This degeneracy renders the measurement unstable if these axes are chosen. However, the minor axis is nondegnerate and can be reliably used to measure the tilt angle of the DM halo with respect to the disk plane.  \nIn 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 DM halo within 50 kpc. Following the GSE merger, the simulation halo appears an oblate profile, in alignment with Iorio & Belokurov (2019), and maintains a negative angular momentum, indicating a retrograde halo with respect to the disk. Joshi et al. (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. For a simple heuristic argument that dynamical friction tries to bring the angular momentum vectors of the halo and disk into alignment, which means the warp will be strongly damped in co-rotating halo but enhenced in counter-rotating halo (Nelson & Tremaine 1995). Thus, the retrograde halo resulting from GSE merger could sustain a long-lived warp.  \nThe 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 Han et al. (2023), they fixed the halo's tilt angle at 25 · and found that the warp need 1.5 Gyr to reach a steady-state amplitude 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 Han et al. (2023), should result in an increasing warp amplitude. However, our simulation could produce warp with a low amplitude at t = 5 . 2 Gyr and show 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.  \nSuch 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  \nthe body-fixed principal axes for the DM halo, resulting in a moment of inertia tensor in the simplified form of [ I ] = diag ( I xx , I yy , I zz ). Consequently, the angular momentum vector can be expressed as H G = { 0 , 0 , H G } .  \nTo 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 is then characterized by changes in these three Euler angles, which can be calculated as follows (Peraire & Widnall 2009):  \n˙ ϕ = H G ( cos 2 ψ I yy + sin 2 ψ I xx ) . (7)  \n˙ θ = H G ( 1 I xx -1 I yy ) sinθsinψcosψ. (8)  \n˙ ψ = H G ( 1 I zz -cos 2 ψ I yy -sin 2 ψ I xx ) cosθ. (9)  \nThe ψ, ϕ, θ are described as 'roll', 'yaw' and 'pitch' angle, and ˙ ψ, ˙ ϕ, ˙ θ are spin, precess and nutate rate respectively. Under the condition that the DM halo exhibits an oblate profile, we align the minor axis with the z-axis in the bodyfixed coordinate system, the evolution of three velocities are shown in Figure 6. The precession rate of DM halo is along the original Z axis. Given that the angular momentum of the dark matter (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 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 single major merger model could not simultaneously fit the both the GSE and the disk warp very accurately.  \nIn 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 < ∼ 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 Milky Way, it's 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.  \nIn Section 5, we find our disk is embedded in a time-evolving tilted and retrograde DM halo. In 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 Gyrs to reach 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 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, lopside, exhibiting both prograde and retrograde precession rate  \nFigure 7. Top left panel shows the orbital decay profile, the merger completed within the first 3 Gyrs. Top middle panel shows all-sky density profile of GSE (blue points/lines), there are two breaks in our simulation and the red dash line is observational stellar halo profile (Deason et al. 2014). Top right panel is the anisotropy parameter β of the GSE debris with Galactocentric distance in a spherical coordinate system, comparing with observational data. Bottom panels are the most representative characteristic of GSE potted with contour profile, the 'sausage' structure is in the left panel and the distribution of the particle binding energy vs. angular momentum in the right panel. Here the negative value of V ϕ 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.  \n<!-- image -->  \nafter 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 appears to sustain the warp.  \nThis 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 the 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 the Yunnan Univ. in China .", 'A. CHARACTERISTIC OF SIMULATED GSE': "Here we list a summary of our GSE model in Figure 7. While most of the galaxy model parameters are adapted from 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, the power-law coefficients are given by ρ ∝ r α gal profile. As 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  \nbreaks in r gal = 15kpc and r gal = 30kpc respectively. The break at 15 kpc is a critical feature of the GSE merger (Naidu et al. 2021) and it also roughly the radius where the disk warp reaches maximum amplitude, which may hint 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 (Deason et al. 2014), and the third coefficients after the second break fits well with that of Ye et al. (2023), as they report that the second break at 27 . 18 kpc with α = 2 . 86. Meanwhile, the top-right panel displays the anisotropy parameter β of GSE debris with Galactocentric distance in a spherical coordinate system, which defined as:  \nβ = 1 -σ 2 θ + σ 2 ϕ 2 σ 2 r . (A1)  \nwhere β quantifies the degree of velocity anisotropy of a system of stellar orbits, where σ i are the velocity dispersions in spherical coordinates. The observational data is from Wu et al. (2022), derived from the LAMOST (Zhao et al. 2006, 2012; Luo et al. 2012; Cui 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 relatively stable to decreasing trend between approximately 17 -28 kpc, and increasing after that. The bottom panels show most prominent features of GSE in both V ϕ -V r panel and the E -L Z panel. We can find the 'sausage' structure in velocity panel which was first discovered in Belokurov et al. (2018) as 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 also reflected in E -L Z panel, where the majority of stars are situated in the central region, and the profile similarly trends toward negative values.", 'REFERENCES': "Abdurro'uf, Accetta K., Aerts C. et al., 2022, ApJS, 259, 35 Amˆores, E. B., Robin, A. C., & Reyl'e, C. 2017, A&A, 602, A67  \nBailin, J. 2003, ApJL, 583, L79  \nBarnes J. E., 2002, MNRAS, 333, 481  \nBattaner, E., & Jim'enez-Vicente, J. 1998, A&A, 332, 809  \nBelokurov, V., Erkal, D., Evans, N. W., Koposov, S. E.,  \nDeason, A. J. 2018, MNRAS, 478, 611  \nBelokurov V., Vasiliev E., Deason A. J., Koposov S. E.,  \nFattahi A., Dillamore A. M., Davies E. Y., Grand R. J. J., 2023, MNRAS, 518, 6200  \nBehroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143  \nBinney J., Jiang I.-G., Dutta S., 1998, MNRAS, 297, 1237 Bignone, L. A., Helmi, A., Tissera, P. B. 2019, ApJL, 883, L5  \nBland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 Bosma, A. 1991, Warped and Flaring HI Disks (Cambridge: Cambridge Univ. Press), 181  \nCabrera-Gadea, M., Mateu, C., Ramos, P., et al. 2024, MNRAS, 528, 4409  \nCiuc˘a I., Kawata, D., Ting, Y.-S., et al. 2024,  \nMNRAS.528L.122C  \nCui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, RAA, 12, 1197  \n- Chen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320 Cheng, X., Anguiano, B., Majewski, S. R., et al. 2020, ApJ, 905, 49 Chrob'akov'a, ˇ Z ., & L'opez-Corredoira, M. 2021, ApJ, 912, 130 Coe D., 2010, preprint (arXiv:1005.0411) Conroy, C., Bonaca, A., Cargile, P., et al. 2019a, ApJ, 883, 107 Cox T. J., Jonsson P., Somerville R. S., Primack J. R., Dekel A., 2008, MNRAS, 384, 386 Deason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30 Dehnen, W., Semczuk, M., & Schonrich, R. 2023, MNRAS, 523, 1556 Dubinski J., & Chakrabarty D., 2009, AJ. 703, 2068-2081 Einasto, J. 1965, Trudy Inst. Astroz. Alma-Ata, 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  \nGaia Collaboration et al. Gaia Data Release 2. Mapping the Milky Way disc kinematics. Astron. Astrophys. 616, A11.  \nWang H., Lopez-Corredoira /acute.ts1 M., Carlin J. L., Deng L., 2018, MNRAS, 477, 2858 Wang J., Hammer F., Athanassoula E., Puech M., Yang Y., Flores H., 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  \nWu, W.B., Zhao, G., Xue,X.X., et al. 2022, ApJ, 924, 23W  \nYe D., Du C., Shi J., Ma J., 2023, MNRAS, 525, 2472 Zhao, G., Chen, Y.-Q., Shi, J.-R., et al. 2006, ChJAA, 6,  \n265 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. doi:10.48550/arXiv.2402.15782"}