A Potential Dynamical Origin of The Galactic Disk Warp: The Gaia-Sausage-Enceladus Major Merger
Mingji Deng Cuihua Du Yanbin Yang Jiwei Liao Dashuang Ye
























































0000-0002-0592-7660]Mingji Deng
School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, P.R. China

:0000-0002-3954-617X]Cuihua Du
School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, P.R. China

0000-0001-7949-3407]Yanbin Yang
GEPI, Observatoire de Paris, Universite PSL, CNRS, Place Jules Janssen, 92195 Meudon, France

0000-0001-6762-5599]Jiwei Liao
School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, P.R. China


School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, Beijing 100049, P.R. China

Cuihua Du

ducuihua@ucas.ac.cn














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




























§ INTRODUCTION
 

   The most common asymmetrical structure in many disk galaxies is the disk warp, as documented in previous studies <cit.>. The Milky Way, a typical example of disk and spiral galaxies, also exhibits a clear disk warp, confirmed by various works <cit.>. Many mechanisms of the Galactic warp have been proposed, including: inflow of intergalactic matter into the halo <cit.>; inflow directly onto the Galactic disk <cit.>; magnetic fields that exist between galaxies <cit.>; interaction between satellite galaxies like Sagittarius <cit.> or Magellanic Clouds <cit.> with the disk; and the bending instability and self-excited warps or internally driven warps of the Galactic disk <cit.>. 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 <cit.>. 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.

   Lines 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 <cit.>. Consequently, the GSE has been studied qualitatively via analogs in cosmological simulations <cit.>, MW zoom-ins <cit.>, and existing merger simulations <cit.>. These studies have achieved success in demonstrating how this merger produces the eccentric debris and reshapes the early Milky Way disk. 
   
   Building upon these findings, <cit.> made a tailored model for this merger by producing a grid of 500 idealized galaxy merger simulations with the GADGET code <cit.> to identify a fiducial model that best matches the H3 survey data <cit.>. 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 <cit.> in their analysis of the age-metallicity relationship from the APOGEE-2 DR17 survey <cit.>. 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 <cit.>, 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. 
   
   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 <cit.>. However, these studies generally prefer a recent interactions with satellite galaxies that drive transient warps. <cit.> 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 <cit.>, which means the warp might be formed in a universal way by galaxies merger.
    
   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 <ref>, we provide a description of the simulation model and the initial conditions. In Section <ref>, we illustrate the long-lived and nonsteady warp model's evolution over time, comparing it with observational data at present. In Section <ref>, we calculate the warp precession rate and demonstrate its lopsided feature. In Section <ref>, 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 <ref> and Section <ref>.



§ MODEL AND SIMULATION
 

   



   In <cit.>, 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 observed 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 <cit.>. 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 <cit.>, 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 <cit.>. 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.

   The numerical initial conditions (ICs) are generated from the disk initial conditions environment (DICE) code <cit.>. 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 <cit.>. 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 <cit.> to replace the Hernquist profile <cit.> used in <cit.>. The Einasto profile, which is similar to the Sérsic profile <cit.> but used for 3-D mass density <cit.>, 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 <ref>, 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 <cit.>. We set the gas disk scale length for the Milky Way and GSE progenitors to 3R_s and 4R_s (R_s is the scale radius of the stellar disk), respectively. These choices follow model assumptions made in <cit.> and <cit.> 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 <cit.>, in which we have implemented the module of star formation and feedback processes as in <cit.>, 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 <cit.>). 

   As these structural choices play a secondary role to the total mass and orbital parameters described in <cit.>, they modeled the orbit with circularity: η=0.5, inclination: θ=15^∘. In <cit.>, they conducted a similar simulation broadly following <cit.>, but changed the inclination to θ=30^∘. Based on the fiducial model of <cit.>, 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 <cit.>, 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 <cit.> 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 <cit.>. It will take several Gyrs for subsequent virialization phase to rebuild the disk <cit.>, 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 <ref>.



§ 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 <ref>. As the GSE merger is estimated to have occurred approximately about 8-11 Gyrs ago <cit.>. In Figure <ref>, 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 <cit.>, 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 <cit.>. The geometric shape of the warp can be approximated by a power-law warp mode:
      
    Z_w(R≥ R_w) =a(R-R_w)^b sin(ϕ -ϕ_w).

      
    Z_w(R< R_w) =0.

   
   Here, 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 <cit.> 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 <cit.> 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. 
   
   In the upper panel of Figure <ref>, 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 <cit.>. 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 H1 data from <cit.> 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 <cit.>.

   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  <cit.>. Accordingly, we plot horizontal lines to represent the height of |Z|=2.5 kpc in the bottom panel of Figure <ref>. 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.
          
   Figure <ref> illustrates the time evolution of the stellar/gas disks warp's amplitude, calculated by fitting Eq. <ref> to the warp every 0.01 Gyr since t=3 Gyr that 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 fitting with stars younger than 1 Gyr. Both disks exhibit similar evolution in Figure <ref>, 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 <cit.>, 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 <ref>, indicating that the warp is a long-lived and nonsteady structure.





§ KINEMATIC WARP MODEL
 



   


   In section <ref>, 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 <cit.>, described as:
      
    V_Z(R,ϕ,t=0)=(V_ϕ/R -ω(R))h_w(R)cos(ϕ-ϕ_ω)+∂ h_w/∂ tsin(ϕ-ϕ_w).

      
   Where h_w(R)=a(R-R_w)^b, as described in Eq. <ref>, 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 <cit.>:
      
    V_Z(R,ϕ)=(V_ϕ/R -ω)h_w(R)cos(ϕ-ϕ_w).

      
   With this kinematic model, after determining the geometric shape through fitting with Eq. <ref>, we can plug in the parameters of h_w and ϕ into Eq. <ref> and measure the procession rate ω using kinematic data from the simulation. Here we select the young stars (age<1Gyr) 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 <ref>. 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∼5km s^-1 kpc^-1. A significant shift in the precession rate is observed around t=5.5 Gyr in Figure <ref>, which can be attributed to the exceptionally low warp amplitude during this period. Comparing both the evolution trend of precession and amplitude in Figure <ref>, we could find that both two consistently reach their extremum simultaneously, indicating a synchronous evolution. In the middle panel of Figure <ref>, 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. <ref>, 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 <cit.>. In the bottom panel of Figure <ref>, refer to <cit.> 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 <cit.>. 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.
   
   In <cit.>, precession was calculated using 12 million red giant stars from Gaia DR2 <cit.>, yielding a precession rate of 10.86± 0.03_stat±3.20_syst km s^-1 kpc^-1 in the direction of Galactic rotation. Recently, the similar model was employed in <cit.>, using 134 Cepheids and combined with the line-of-sight velocity (RV) from Gaia DR3 <cit.>, determining a lower warp precession rate of 4.9±1.6 km s^-1 kpc^-1, our simulation is very close to these observation in different ranges. Another study by <cit.> 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 <ref>, 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 <cit.> for R<∼15 kpc. Both <cit.> and <cit.> 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 <cit.>, which might make precession rate rise in observation beyond 15 kpc.

   
   


   Theoretical 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 oblate-like halo should produce a retrograde warp precession <cit.>. However, warps in oblate halos are not appropriate to explain the observed warps or they could not sustain a prominent warp <cit.>. Typical precession periods from such models or those produced by a misaligned outer torus of later accreted material <cit.>, are approximately between 4 Gyr and 40 Gyr, corresponding to warp precession rates between 1.5 km s^-1 kpc^-1 and 0.1 km s^-1 kpc^-1. In Figure <ref>, the precession rate determined after the merger significantly surpasses all predicted values throughout the simulation 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 <cit.>, 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.





§ TILTED AND RETROGRADE DM HALO
 


   


     
   


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

    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 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 <cit.>:
    
    r_i^2∝-I_i+I_j+I_k/2.

    Once we obtain the moments of inertia, we can calculate the major-to-intermediate and major-to-minor axes ratios as follows:
    
    r_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).
  
    
    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 <cit.> and <cit.>. 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.

    In the middle and bottom panels of Figure <ref>, 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 <cit.>, and maintains a negative angular momentum, indicating a retrograde halo with respect to the disk. <cit.> 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 <cit.>. Thus, the retrograde halo resulting from GSE merger could sustain a long-lived warp. 
    
    The top panel of Figure <ref> 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 <cit.>, 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 <cit.>, 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 <ref>. 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]=diag(I_xx,I_yy,I_zz). Consequently, the angular momentum vector can be expressed as  𝐇_𝐆={  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_ZR_XR_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 <cit.>:
    
    ϕ̇ = H_G(cos^2ψ/I_yy+ sin^2ψ/I_xx).

    
    θ̇ = H_G(1/I_xx - 1/I_yy)sinθ sinψ cosψ.

    
    ψ̇ = H_G(1/I_zz - cos^2ψ/I_yy-sin^2ψ/I_xx )cosθ.



     The ψ,ϕ,θ 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 body-fixed coordinate system, the evolution of three velocities are shown in Figure <ref>. 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 <ref>, where the precession remains predominantly retrograde throughout most of the period.
    




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

    In Section <ref>, 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.

    In Section <ref>, we find our disk is embedded in a time-evolving tilted and retrograde DM halo. In <cit.>, 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.




§ 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 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 appears to sustain the warp. 




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























§ CHARACTERISTIC OF SIMULATED GSE
 
      
       
    

    Here we list a summary of our GSE model in Figure <ref>. While most of the galaxy model parameters are adapted from <cit.>, 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 <cit.> 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 <cit.> 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 proﬁle between two breaks is a good match to that found for the inner halo <cit.>, and the third coefﬁcients after the second break fits well with that of <cit.>, 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:
      
    β = 1-σ_θ^2+σ_ϕ^2/2σ_r^2.

    where β quantiﬁes 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 <cit.>,  derived from the LAMOST <cit.> 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 <cit.> 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.
    
    












  [Abdurro’uf et al.(2022)]AbdurroufAbdurro’uf, Accetta K., Aerts C. et al., 2022, ApJS, 259, 35
  [Amôres et al.(2017)]AmoresAmôres, E. B., Robin, A. C., & Reylé, C. 2017, A&A, 602, A67
  [Bailin(2003)]BailinBailin, J. 2003, ApJL, 583, L79
  [Barnes(2002)]BarnesBarnes J. E., 2002, MNRAS, 333, 481
  [Battaner & Jiménez-Vicente (1998)]BJBattaner, E., & Jiménez-Vicente, J. 1998, A&A, 332, 809
  [Belokurov et al.(2018)]Belokurov Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., Deason, A. J. 2018, MNRAS, 478, 611
  [Belokurov et al.(2023)]Belokurov23Belokurov V., Vasiliev E., Deason A. J., Koposov S. E., Fattahi A., Dillamore A. M., Davies E. Y., Grand R. J. J., 2023, MNRAS, 518, 6200
  [Behroozi et al.(2019)]BehrooziBehroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  [Binney et al.(1998)]BinneyBinney J., Jiang I.-G., Dutta S., 1998, MNRAS, 297, 1237

  [Bignone et al.(2019)]BignoneBignone, L. A., Helmi, A., Tissera, P. B. 2019, ApJL, 883, L5
  [Bland-Hawthorn & Gerhard(2016)]BHGBland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
  [Bosma(1991)]BosmaBosma, A. 1991, Warped and Flaring HI Disks 
  (Cambridge: Cambridge Univ. Press), 181
  [Cabrera-Gadea et al.(2024)]CCCabrera-Gadea, M., Mateu, C., Ramos, P., et al. 2024, MNRAS, 528, 4409
  [Ciucă et al.(2024)]CiucaCiucă I., Kawata, D., Ting, Y.-S., et al. 2024, MNRAS.528L.122C
  [Cui et al.(2012)]CuiCui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, RAA, 12, 1197
  [Chen et al.(2019)]ChenChen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320
  [Cheng et al.(2020)]ChengCheng, X., Anguiano, B., Majewski, S. R., et al. 2020, ApJ, 905, 49
  [Chrobáková et al.(2021)]ChrChrobáková, Ž., & López-Corredoira, M. 2021, ApJ, 912, 130
  [Coe(2010)]CoeCoe D., 2010, preprint (arXiv:1005.0411)
  [Conroy et al.(2019a)]ConroyConroy, C., Bonaca, A., Cargile, P., et al. 2019a, ApJ, 883, 107
  [Cox et al.(2008)]CoxCox T. J., Jonsson P., Somerville R. S., Primack J. R., Dekel A., 2008, MNRAS, 384, 386
  [Deason et al.(2014)]DeasonDeason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30
  [Dehnen et al.(2023)]Dehnen23Dehnen, W., Semczuk, M., & Schönrich, R. 2023, MNRAS, 523, 1556
  [Dubinski, J. & Chakrabarty(2009)]DubinskiDubinski J., & Chakrabarty D., 2009, AJ. 703, 2068–2081
  [Einasto(1965)]EinastoEinasto, J. 1965, Trudy Inst. Astroz. Alma-Ata, 51, 87 
  [Elias et al.(2020)]EliasElias, L. M., Sales, L. V., Helmi, A., Hernquist, L. 2020, MNRAS, 495, 29
  

  [Fattahi et al.(2019)]FattahiFattahi, A., Belokurov, V., Deason, A. J., et al. 2019, MNRAS, 484, 4471
  [Freudenreich et al.(1994)]FreudenreichFreudenreich, H. T., Berriman, G. B., Dwek, E., et al. 1994, ApJL, 429, L69
  
  [Gaia Collaboration et al.(2018)]GaiaGaia Collaboration et al. Gaia Data Release 2. Mapping the Milky Way disc kinematics. Astron. Astrophys. 616, A11.
  [Gaia Collaboration et al.(2023a)]Gaia3Gaia Collaboration, Drimmel, R., Romero-G´omez, M., et al. 2023a, A&A, 674, A37
  [Grand et al.(2017)]Grand17Grand R. J. J., Gómez F.,  Marinacci F., et al. 2017, MNRAS, 467, 179
  [Grand et al.(2020)]GrandGrand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603
  [Hammer et al.(2009)]HammerHammer F., Flores H., Puech M., Yang Y. B., Athanassoula E., Rodrigues M., Delgado R., 2009, A&A, 507, 1313
  [Han et al.(2023)]Han23Han, J. J., Conroy, C., Hernquist, L., 2023NatAs...7.1481H
  [Han et al.(2022)]Han22Han, J. J., Conroy, C., Johnson, B. D., et al. 2022a, AJ, 164, 249
  [He(2023)]HeHe. 2023, ApJ, 954L, 9H
  [Helmi et al.(2018)]HelmiHelmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85
  [Hernquist(1990)]HernquistHernquist, L. 1990, ApJ, 356, 359
  
  [Hopkins et al.(2009)]Hopkins09Hopkins P. F., 2009, in Jogee S., Marinova I., Hao L., Blanc G. A., eds, ASP Conf. Ser. Vol. 419, Galaxy Evolution: Emerging Insights and Future Challenges. Astron. Soc. Pac., San Francisco, p. 228
  [Hopkins(2015)]Hopkins15Hopkins P. F., 2015, MNRAS, 450, 53
  [Hopkins et al.(2018)]Hopkins18Hopkins P. F., Wetzel A.,  Kereš D., 2018b, MNRAS, 480, 800
  [Huang et al.(2024)]HuangHuang, Y. et al. A slightly oblate dark matter halo revealed by a retrograde precessing galactic disk warp. Nature Astronomy early access (2024).
  [Ideta et al.(2000)]IdetaIdeta M., Hozumi S., Tsuchiya T., Takizawa M., 2000, MNRAS, 311, 733
  [Iorio & Belokurov(2019)]IorioIorio, G. & Belokurov, V. 2019, MNRAS, 482, 3868
  [Jeon et al.(2009)]JeonJeon, M., Kim, S. S., & Ann, H. B. 2009, ApJ, 696, 1899
  [Jiang & Binney(1999)]JBJiang, I.-G., & Binney, J. 1999, MNRAS, 303, L7
  [Joshi et al.(2024)]JoshiJoshi, R., Widrow, Lawrence M. 2024MNRAS.527.7781J
  [Kazantzidis et al.(2004)]KazantzidisKazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJL, 611, L73
  [Kerr(1957)]KerrKerr, F. J. 1957, AJ, 62, 93

  [Koppelman et al.(2020)]KoppelmanKoppelman, H. H., Bos, R. O. Y., Helmi, A. 2020, A&A, 642, L18
  [Laporte et al.(2019)]LaporteLaporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134
  [Lemasle et al.(2022)]LemasleLemasle et al. 2022, A&A, 668, A40
  [Levine et al.(2006)]LevineLevine E. S., Blitz L., Heiles C., 2006, Science, 312, 1773
  [Luo et al.(2012)]LuoLuo, A. L., Zhang, H.-T., Zhao, Y.-H., et al. 2012, RAA, 12, 1243
  [López-Corredoira et al.(2002a)]LCLópez-Corredoira, M., Betancort-Rijo, J., & Beckman, J. E. 2002a, A&A, 386, 169
  [López-Corredoira et al.(2014)]LC14López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101
  [Li et al.(2023)]Li23Li, X., Wang, H.-F., Luo, Y.-P., et al. 2023, ApJ, 943, 88
  [Metropolis et al.(1953)]MetropolisMetropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087
  

  [Naidu et al.(2020)]Naidu2020Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48
  [Naidu et al.(2021)]Naidu2022Naidu R. P., Conroy C., Bonaca, A., et al. 2021, ApJ, 923, 92
  [Nelson & Tremaine(1995)]ntNelson R. W., Tremaine S., 1995, MNRAS, 275, 897
  [Ostriker & Binney(1989)]OBOstriker, E. C., & Binney, J. J. 1989, MNRAS, 237, 785
  [Peraire & Widnall(2009)]Peraire Peraire J., Widnall S., 16.07 Dynamics, Fall 2009 Version 2.0
  [Perret et al.(2014)]Perret14Perret, V., Renaud, F., Epinat, B., et al. 2014, A&A, 562, A1
  [Perret(2016)]Perret16Perret, V. 2016, DICE: Disk Initial Conditions Environment, Astrophysics Source Code Library, ascl:1607.002
  [Poggio et al.(2020)]PoggioPoggio, E., Drimmel, R., Andrae, R., et al. 2020, Nature Astronomy, 4, 590
  [Quinn & Binney(1992)]QBQuinn, T., & Binney, J. 1992, MNRAS, 255, 729

  [Revaz & Pfenniger(2004)]RPRevaz, Y., & Pfenniger, D. 2004, A&A, 425, 67

  [Reshetnikov & Combes(1998)]ReshetnikovReshetnikov, V., & Combes, F. 1998, A&A, 337, 9
  [Reylé et al.(2009)]ReyleReylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819
  [Rocha et al.(2008)]RochaRocha M., Jonsson P., Primack J. R., Cox T. J., 2008, MNRAS, 383, 1281
  [Rodrigues et al.(2012)]RodriguesRodrigues M., Puech M., Hammer F., Rothberg B., Flores H., 2012, MNRAS, 421, 2888
  [Roskar et al.(2010)]RoskarRoskar, R., Debattista, V. P., Brooks, A. M., et al. 2010, MNRAS, 408, 783
  [Sánchez-Saavedra et al.(1990)]SS90Sánchez-Saavedra, M. L., Battaner, E., & Florido, E. 1990, MNRAS, 246, 458
  [Sánchez-Saavedra et al.(2003)]SS03Sánchez-Saavedra, M. L., Battaner, E., Guijarro, A., López-Corredoira, M., & Castro-Rodríguez, N. 2003, A&A, 399, 457
  [Sellwood & Debattista(2022)]SDSellwood, J. A., & Debattista, V. P. 2022, MNRAS, 510, 1375
  [Sérsic(1968)]SersicSérsic, J. L. 1968, Atlas de galaxias australes (Cordoba, Argentina: Observatorio Astronomico, 1968)
  [Shao et al.(2021)]ShaoShao, S., Cautun, M., Deason, A., & Frenk, C. S. 2021, MNRAS, 504, 6033
  [Shen & Sellwood(2006)]ShenShen, J. & Sellwood, J. A., 2006, MNRAS, 370, 2
  [Skowron et al.(2019a)]SkowronSkowron, D. M., Skowron, J., Mróz, P., et al. 2019a, AcA, 69, 305
  [Snaith et al.(2014)]SnaithSnaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJ, 781, L31
  [Springel et al.(2005)]Springel05Springel, V. 2005, MNRAS, 364, 1105
  [Springel et al.(2021)]Springel21Springel, V., Pakmor, R., Zier, O., Reinecke, M. 2021, MNRAS, 506,2871
  [Stelea et al.(2024)]SteleaStelea, I. A., Hunt, J. A. S., & Johnston, K. V. 2024, The Milky Way’s rowdy neighbours: The effects of the Large Magellanic Cloud and Sagittarius Dwarf on the Milky Way Disc
  [Wang et al.(2018)]WangWang H., Lopez-Corredoira ´ M., Carlin J. L., Deng L., 2018, MNRAS, 477, 2858
  [Wang et al.(2012)]WangjWang J., Hammer F., Athanassoula E., Puech M., Yang Y., Flores H., 2012, A&A, 538, A121
  [Weinberg & Blitz(2006)]WBWeinberg, M. D., & Blitz, L. 2006, ApJL, 641, L33
  [White & Frenk(1991)]WhiteWhite, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52
  [Wu et al.(2022)]WuWu, W.B., Zhao, G., Xue,X.X., et al. 2022, ApJ, 924, 23W
  [Ye et al.(2023)]YeYe D., Du C., Shi J., Ma J., 2023, MNRAS, 525, 2472
  [Zhao et al.(2006)]Zhao06Zhao, G., Chen, Y.-Q., Shi, J.-R., et al. 2006, ChJAA, 6, 265
  [Zhao et al.(2012)]Zhao12Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, RAA, 12, 723
  [Zhou et al.(2024)]ZhouZhou, X.Y., Chen, X.D., Deng, L.C., et al. doi:10.48550/arXiv.2402.15782












