Mapping the Galactic disk with the LAMOST and Gaia Red clump sample. VI: An evidence for the long lived non-steady warp of non-gravitational scenarios

By combining LAMOST DR4 and Gaia DR2 common red clump stars with age and proper motion, we analyze the amplitude evolution of the stellar warp independently of any assumption with a simple model. The greatest height of the warp disk increases with Galactocentric distance in different populations and it is dependent on the age: the younger stellar populations exhibit stronger warp features than the old ones, accompanied with the warp amplitude $\gamma ({\rm age})$ decreasing with age and its first derivative $\dot{\gamma} ({\rm age})$ is different from zero. The azimuth of line of nodes $\phi _w$ is stable at $-$5 degree without clear time evolution, which perfectly confirms some of previous works. All these self-consistent evidences support that our Galactic warp should most likely be a long lived, but non-steady structure and not a transient one, which is supporting the warp is originated from gas infall onto the disk or other hypotheses that suppose that the warp mainly affects to the gas, and consequently younger populations tracing the gas are stronger than the older ones. In other words, the Galactic warp is induced by the non-gravitational interaction over the disk models.


INTRODUCTION
Most spiral galaxies have a warped disk (Reshetnikov et al. 1998;Sánchez-Saavedra et al. 1990, 2003, although in some spirals warps cannot be observed due to their low inclination. Like many spiral galaxies in the universe, Milky Way's (MW) warp was detected by neutral hydrogen (HI) gas many years ago (Kerr 1957;Bosma 1981;Briggs et al. 1990;Levine et al. 2006b;Nakanishi et al. 2003). Dust has also been observed in Marshall et al. (2006). Furthermore, there are several works that measure parameters of the stellar Galactic warp from density maps/star counts (López-Corredoira et al. 2002b;Momany et al. 2006;Reylé et al. 2009;Amôres et al. 2017). There is also an intuitive 3 dimensional (3D) map of the Galactic warp traced by classical Cepheids (Chen et al. 2019;, including precession measurements consistent with the Brigg's rule (Briggs et al. 1990) and showing us that its amplitude of northern part is very prominent and stronger than that of the southern part (Skowron et al. 2019b). Kinematic signatures of the Galactic warp were studied by Huang et al. (2018); Schönrich et al. (2017) in the vertical velocity, angular momentum, azimuthal velocity or guiding center radius parameters space, showing that the Galactic warp is not a transient struc-ture, which is consistent with the simple calculation of Wang et al. (2020a). More kinematical signals could also be found in Smart et al. (1998); Poggio et al. (2018); Romero− Gómez et al. (2019), etc. The clear differences between thin disk and thick disk warp classified by metallicity and abundance is shown in one of our series of works (Li et al. 2020).
The mechanisms of formation of the gas and stellar warp were proposed in many works. Debattista & Sellwood (1999); Shen et al. (2006) showed that warps were produced by the dynamical friction between a misaligned rotating halo and disk. In some cases, this misalignment might be related to misaligned gas infall (Ostriker et al. 1989;Quinn et al. 1992;Bailin et al. 2003b). Burke et al. (1957); Weinberg et al. (2006) suggested it was caused by the interaction with the Magellanic Clouds and Bailin et al. (2003a) proposed the cause of interaction with Sagittarius (Bailin et al. 2003a) ;Hunter &Toomre (1969) claimed that Magellanic Clouds mass is not enough to explain the observed amplitude of the warp, however, new elements of amplification have been introduced by Weinberg et al. (2006). Other scenarios, that is to say, perturbations by dwarf satellites (Weinberg 1995;Shen et al. 2006), intergalactic magnetic fields influence (Battaner & Jiménez-Vicente 1998) or accretion/infall of the intergalactic medium flows directly onto the disk (López-Corredoira et al. 2002a), disk bending instabilities by Revaz et al. (2004), etc., also appeared. These and other mechanisms were proposed but some observational evidences may favour some of them, the different interpretations on the formation of the warp are still being hotly debated. In any case, we know the kinematical distributions of vertical bulk motions will be contributed by warp asymmetrical structure (Wang et al. 2018a(Wang et al. , 2020a, in return, vertical motions can be used to constrain the warp's properties. Warp in the Milky Way bends up upwards and downwards in the north and south hemisphere separately with different amplitude at least in the gas . The amplitude of the warp clearly increases strongly with radius and varies with the azimuth angle (Li et al. 2020;Liu et al. 2017c;López-Corredoira et al. 2014). A linear simple relation between the amplitude of the warp and the Galactocentric distance was used to explain the increase trend of vertical velocity with vertical angular momentum (Schönrich et al. 2017). In López-Corredoira et al. (2014), the vertical bulk motions are contributed by a warp with modeling as a set of circular rings that are rotated and whose orbit is in a plane with angle with respect to the Galactic plane, with vertical amplitude z w (R, φ) = γ(R-R ) α sin(φ-φ w ), where R is the Galactocentric distance and φ is the Galactocentric azimuth. With the calculation of γ,γ by assuming φ w and α are constant, their work indicated that most likely the main S-shaped structure of the warp is a long-lived feature.
As mentioned in López-Corredoira et al. (2014), the precision of vertical velocity can be increased by at least an order of magnitude with the help of the Gaia proper motion (Gaia Collaboration et al. 2018;Gaia Collaboration: Katz et al. 2018) together with spectroscopically classified red clump stars, e.g. LAMOST (Cui et al. 2012;Deng et al. 2012;Liu et al. 2014;Zhao et al. 2012). In addition to the unprecedented proper motion, we have stars age nowadays, so we have this first chance to research the evolution of the warp structure properties. For this work, we are motivated to use López-Corredoira et al. (2014) simple model and LAM-OST (Large Sky Area Multi-Object Fibre Spectroscopic Telescope) DR4, and Gaia DR2 common red clump stars to investigate the warp amplitude, its first derivative and greatest height variation with different age populations so that we could get better constraint for the γ ,γ, and φ w . Hence, we could offer our interpretation of the formation and evolution history of Milky Way warp. This paper is structured as follows: The section 2 is about how we select our red clump stars, velocity derivation, and the vertical velocity distribution in different age populations. Model and method are introduced in the section 3 . Our results will be shown in section 4 and discussions are displayed in the section 5. Finally we will give our conclusions of this work.

THE SAMPLE SELECTION
During this work, we use the red clump giants selected from the LAMOST Galactic spectroscopic surveys and Gaia astrometric survey. The scientific motivations and target selections of LAMOST phase I can be found in Cui et al. (2012); Deng et al. (2012);Liu et al. (2014); Zhao et al. (2012). Now, we are entering into the Phase-II. Its fiber is 3.3 arcsec and the mean seeing during LAMOST observations is around 3 arcsec and the spatial resolution of LAMOST should be around 5 arcsec. Selection functions of LAMOST is almost a flat along apparent magnitude, more details can also be found in Carlin et al. (2012); Yuan et al. (2015); Liu et al. (2017a).
The red clump stars selection details can be found in Huang et al. (2020Huang et al. ( , 2015, the distance and age is determined by the Kernel Principal Component Analysis (KPCA) method, which could be found more details in Xiang et al. (2017b,c). According to Huang et al. (2020), the distance uncertainties are around 5-10% and age uncertainties are around 30%, which have been used quite well in Wang et al. (2019Wang et al. ( , 2020a, and Red clump stars are well known horizontal stars and standard candles so that it is not strange that the distance error of our sample is small. The Gaia DR2 catalogue contains highprecision positions, parallaxes, and proper motions for 1.3 billion sources as well as line-of-sight velocities for 7.2 million stars. For stars of G < 14 mag, the median uncertainty is 0.03 mas for the parallax and 0.07 mas yr −1 for the proper motions. In order to get reliable stellar parameters and try to reduce the halo contamination from the 0.14 million sample mainly consisted of primary red clumps (RCG) with few contaminations of secondary red clumps and Red Giant Branch(RGB) stars, we use the latest sample of this catalog to investigate our scientific target carefully according to the following criterions: 1. Sample without parameters such as distance, radial velocity, temperature and surface gravity are removed.
3. Stars with LAMOST spectroscopic SNR > 20 and age less than 15 Gyr are included. We derive the 3D velocities and Cartesian coordinates on the basis of coordinate transformation described in Galpy (Bovy et al. 2015), which is more convenient and direct than those described on by one in Wang et al. (2018aWang et al. ( , 2019Wang et al. ( , 2020a. The vertical angular momentum distributions of the sample associated with the error analysis in the R, Z plane of the Cartesian coordinate system are shown in left panel of Fig. 1, the mass determination is also given in this catalog with the uncertainties of 15%. Meanwhile, the projected age, [F e/H] and its measured error distributions are displayed in the middle and right panel respectively. It reasonably implies that the angular momentum of stars increases with radial distance in the disk, including the corresponding errors, due to that disk total momentum is approximate to be vertical momentum. The age distribution for which the stars inside are relatively older and those outside are generally younger support the inside−out formation of the Milky Way disk, age error is smaller comparing with the age value, albeit probably some possible systematic errors that remain to be ignored, it can still give us the good chance to map the dynamical structures in different age populations. As shown in the bottom right one of Fig 1, with Galactic radial distance increasing, the metallicity has a negative trend and the error is also becoming larger with value of around 0.1−0.15 dex, which is reasonable. Here we want to emphasize that, by using this recent updated sample, some asymmetrical structures such as radial or bulk motions reconstructed here are very similar to our previous series of works about the Galactoseismology (Wang et al. 2020a)  We could see the vertical velocity (V Z ) distribution of our sample with radial distance (R) in different age populations with different stellar ages in Fig 2. As shown in each sub-figure, the vertical velocity increases with radial distance in different mono-age populations from 0 to 10 km s −1 , it is definitely reflecting the warp signals. And the vertical velocities of systematic increase for most age bins are around 6−8 km s −1 from 8 to 14 kpc except the last one, which is similar to Poggio et al. (2018) for the value of 5−6 km s −1 . Although there are some oscillations due to the Poissonian noise, as we will mention again in the next section 4, in order to get more points to show how the vertical velocity along with the distance and ensure we have enough data to do fitting, here we plot the velocity profile with the bins containing at least 20 stars. If we enlarge the bin size, the oscillation will be reduced. It is also clear to see that the younger age bins of top three panels have significantly smaller errors compared with the bottom three panels. Again, it might be caused by Poissonian noise. Moreover, the age accuracy is also becoming worse and worse as the stars are older and older. The youngest population vertical velocity is significantly larger than the oldest one in the top left and bottom right panel respectively, which might imply a warp amplitude difference.

MODEL
With the assumption that this vertical motion is contributed by warp, modeled as a set of circular rings that are rotated and whose orbit is in a plane with angle i w (R) with respect to the Galactic plane, as we can see many more details in López-Corredoira et al. (2014), the modeling process is displayed as follows: where φ w is the azimuth of the line of nodes, and z w is the height of the disk over the b = 0 plane. We assume the greatest height of the warp to be and a variable line of nodes which has no extremely slow precession to do fitting, (i.e.φ w γ) and changed with the shape of the warp are adopted. To do so, we also assume a constant rotation speed Ω(R, z) = Ω LSR = 238 km/s; this may be slightly reduced for high R or high |z| (López-Corredoira et al. 2014), but the order of magnitude does not change, and V Z is only weakly affected by a change of the rotation speed. Combining all these formulas and assumptions, we derive, for low angles i w (R), the low height disk warp model can be simplified reasonably as: where φ w is the azimuth of the line of nodes, γ is the amplitude of the warp andγ describes the warp amplitude evolution, that is to say, −d(γ)/d(age). We assume the exponent α = 2 (kpc −1 ) ( During this likelihood, R i is the ith point of the fixed Z grid in different age bins. It is emphasized here that each R i is naturally corresponding to a φ i during the process and warp could vary with radius and azimuth angle, we use the information of Fig 2 to constrain warp for this work. Please notice that we set the parameter asγ = −d(γ)/d(age) by adopting a joint likelihood with 12 parameters to do simultanous fitting of all of age bins. In order to get the convergent parameters and save computer time, we just choose a relatively smaller sampling size in our simulation. For a test, we also set the larger sampling size in MCMC, but the pattern is stable. As an attempt to explore the amplitude, line of nodes and maximum stellar warp height with age, our results are shown in the next section.

Simplified analytical model fitting
It is well known and mentioned that the vertical bulk motions can be excited by the warp (Roškar et al. 2010;López-Corredoira et al. 2014;Wang et al. 2018aWang et al. , 2020a, which implies clearly vertical upward motions can be used to reveal the properties of warp, such as warp amplitude, precession rate and so on. The fitting results of the work in different age bins, by fitting simultaneously all the age bins, with α = 2 in model, are displayed in Fig 3. Some detailed warp features are revealed in the likelihood distribution of the parameters (γ and φ w ) drawn from the MCMC simulation in the next section.

4.2.
Warp parameter γ,γ and φ w evolution with age In Fig. 4, the amplitude γ (kpc −1 ) evolution of warp with age is shown in the top left panel. We use the median value of six age bins, that is to say, [1, 3 ,5, 7, 9, 12] Gyr, as the x axis value; and y axis value is the fitting results calculated by the MCMC according to its Gaussian distribution. The probability distribution and its peak are similar in all of the cases. The error is re-calculated again by bootstrap process. We can see that there is variation of the γ(kpc −1 ) with age with relatively large 1−σ error, all these values are decreasing with age. Correspondingly, the bottom left panel is the warp amplitude derivative variationγ (kpc −1 Gyr −1 ) distribution with age, it has a variable increasing trend, these values are different from zero, implying that the warp is always existing but not a stationary structure (γ (kpc −1 Gyr −1 ) ∼ 0) and there is clear difference of populations existing. Moreover, there is also a stable feature for the azimuth of line of nodes in all populations, the value is almost fixed at about −5 • (degree) for the distribution displayed in the top right panel, the relative larger error makes it to be shown as a flatline.
We are considering the variation of φ w with time negligible. Poggio et al. (2019) assume that φ w is variable while the variation of γ (the amplitude of the warp) is negligible. We cannot distinguish in the vertical motion maps which is the dominant factor in V Z ,γ orφ w , but  the variation of φ w due to precession is too small to be detected in all of the warp models (López-Corredoira et al. 2002a;Dubinski et al. 2009;Jeon et al. 2009). Moreover, the fit depend on the model for warp shape. Poggio et al. (2019) have adopted a too large amplitude warp of young populations, which is not applicable to the old populations of Gaia (Chrobáková et al. 2020), and when they get too high values of the vertical motions without precession (their Fig. 3/bottom/left), they have to introduce a too high precession to compensate it. We are skeptical about the validity of their results based on wrong assumptions. As a natural product, theγ/γ is shown in the bottom right sub-figure, which is different from zero and increases with age. The figures derived from the MCMC simulations for all these Age bins could give us a reliable estimation for these parameters of γ (kpc −1 ) and φ w (deg) thanks to the peak of the maximum likelihood. Note that γ has units of kpc −1 for our adopted value of α=2; it has no units when we set α=1.

Z w : Greatest height of warp for younger ages
In Fig. 5, we show the distribution of the greatest height of warp disk with different age bins (bottom panel) and distance of different populations(top panel). The top panel suggests that there is an increasing trend for the height along with distance in all age bins. We also see the younger populations are higher than the old ones. Note that there are six age bins but two lines are overlapped. For the bottom one, we can clearly see there is clear decreasing pattern for all median heights in all age populations. It is consistent with the results in Fig. 4, meaning that the warp appears to be a long lived not-stationary structure and, more importantly, there is a clear difference for different populations again.

The comparison of Model and Sample
In order to check our robustness of the derived results, we have finished the comparison of Model of Eq. (4) and data in six age populations during this work as shown in Fig. 6. From the top left panel to the bottom right one, the median values of V Z for [1, 3 ,5, 7, 9, 12] Gyr populations are plotted. The black bold solid line with error bars is the observed vertical velocity distribution with Galactocentric radial distance in each panel; the blue one is the model profile with the Monte-Carlo fitting, the cyan one is the model profile plus 1σ and the green is model profile plus −1σ. As we can see here, for the general trend, the matching results of all populations between the data and model are good within the uncertainties. The goodness of fit is not good enough at the left boundary in some populations. However, again, the general trend is quite good in 1σ for most of the points. Few mismatches shown here are caused by some reasons such as, in some regions, our model is a simple model, it is not very perfect and suitable for us to describe different populations in some cases; when the distance is larger, the error of stellar parameter and age precision is becoming worse, making it cannot be fitted quite well; the Poisson noise due to the the number of our sub-sample is not large enough; the Sun is not being on the line-of-nodes thus cause some stars not to move towards anticenter with different directions; extinctions in some regions are not perfect.
In order to try to test the Poissonian noise to mislead us for the conclusion, we just compare our model with data consisted of at least 100 stars in each bin, which has small differences for the pattern shown in Fig. 6 and Fig.  2. The overall trend is better except fewer boundary points deviate, which can not change our conclusion at all. So we suggest all features mentioned in the previous paragraphs are robust and intriguing.
Therefore, in short summary, we think that these stable features are real and these observational evidences strongly support that the warp is a long lived and notsteady S shape structure, and it also implies the warp evolution is relatively uniform and there is a clear difference for different age bins. So far, we only have a few points span from 8−14 kpc with relatively large and not-perfect age error, it would be worth for us to further investigate in more details for this structure with the help of sample consisted of larger distance range, more stars, and more accurate age. Some scenarios are given in the next section.

DISCUSSIONS
Vertical non-axisymmetries and wave-like density patterns are found in the solar neighborhood and in the outer disk (Widrow et al. 2012;Williams et al. 2013;Xu et al. 2015;Wang et al. 2018a,b,c;Carlin et al. 2013;Carrillo et al. 2018Carrillo et al. , 2019Pearl et al. 2017). As mentioned and implied in Wang et al. (2020a), many mechanisms including warp might be coupled together under the same complexed distribution function to cause the complicated Galacto-seismology signals. Scenarios for producing these structures such as warp dynamics, minor mergers or interactions with nearby dwarfs or satellites (Gómez et al. 2013;Laporte et al. 2018;D'Onghia et al. 2016;Laporte et al. 2019) and the effects of even lower-mass dark matter sub-halos have also been invoked as a possible explanation (Widrow et al. 2014). Notice that vertical velocity asymmetry can be applied to constrain the warp signal and it is acceptable that we use the vertical motions to acquire the warp amplitude  Figure 4. Top one shows the amplitude γ (kpc −1 ) evolution with age, the median value of six age bins are adopted in x axis: [1, 3 ,5, 7, 9, 12] Gyr; y axis is the MCMC values. Correspondingly, the bottom left panel is theγ (kpc −1 Gyr −1 ) distribution for age with the error determined by bootstrap error. The top right panel is the stable azimuth of line of node along with age and the corresponding bottom right one is the natural product (γ/γ) based on the left two figures, which is also increasing accompanied with errors, all these figures errors are from bootstrap process. and its variation, although other causes different from warp may contribute to the vertical motions too. The kinematical features of the Galactic warp around its line-of-nodes located close to the Galactic anti-center region are discussed in Liu et al. (2017c), where it was found that the vertical bulk motion of younger red clump stars are significantly larger than that of the older ones, which is consistent with our results. We have got a not-steady warp. A variation of warp amplitude with stellar population age is in principle against a steady warp due to steady gravitational forces and it is more in favour of the models in which gas is necessary for warp formation, or similar scenarios. The reason is that the young population, tracing the gas, will always have larger warp amplitude.
According to Skowron et al. (2019b), there are mainly two classes of warp formation mechanisms. One is that the warp formed by the gravitational interactions, for example, with satellite galaxies or a misaligned dark matter halo. The other one is non-gravitational mechanisms, e.g, accretion of intergalactic gas onto the disc (López-Corredoira et al. 2002a), or interactions with intergalactic magnetic fields. Non-gravitational mechanisms such as magnetic fields or hydrodynamical pres-sure from infalling gas would act on the gas and only affect the young stars (Guijarro et al. 2010;Sellwood et al. 2013), thus we should see all signals of younger ones are stronger than the old ones, and we do get that clearly. Therefore, we think the gravitational scenario should not be the reason at least for this tracer.
Young population traces the gas, whereas old population had more time to reduce the amplitude of the warp due to the self-gravity in the models in which the torque affects mainly to the gas and not to the stars. Our current results support the warp might be contributed by the non-gravitational interaction models, which do not agree with the viewpoints of Poggio et al. (2018) by using upper main sequence stars and giants as two age populations.
An age dependence both on position and kinematics of the Galactic warp has been observed by Amôres et al. (2017); Romero− Gómez et al. (2019), and confirmed for our results in Fig. 2 and Fig. 5. We could see there are some vertical velocity differences and the greatest height difference in different age bins clearly. Romero− Gómez et al. (2019) showed the amplitude of OB stars corresponding to younger populations is weaker than the RGB stars corresponding to older ones by calculating (3) with some assumptions, distributions with Galactic radial distance, in which the different color solid lines represent different ages of [1, 3 ,5, 7, 9, 12] Gyr. It has positive gradient and increased with R in different populations. The bottom one is the median value of the Zw with age, a clear decreasing are here. Notice that there are six age bins but two lines are overlapped each other.
the onset radius and height of warp, and thus suggesting that warping disk of older population is more pronounced or stronger, which is not consistent with our results implying that the warp amplitude is variable and decreasing in different age bins. The discrepancy might be due to the fact that we use the direct age bins and greatest height to describe the warp amplitude, but Romero− Gómez et al. (2019) use the height and indirect age indicators with OB and RGB stars, so the methods, assumptions and population effects will be important for the discrepancy. Furthermore,  showed that their results, by using Cepheids similar with OB stars, are in contradiction to Romero− Gómez et al. (2019) with that the Cepheids warp height is similar with RGB stars. We also notice that the recent results of Chen et al. (2019) about the warp for very young populations give the highest amplitude of the warp so far for stellar populations, which also supports our conclusion that younger ones are stronger ones. A young population warp larger than the old population one has also been demonstrated by Chrobáková et al. (2020) using Gaia-DR2 density maps extended up to R = 20 kpc, thanks to the use of deconvolution techniques of parallax errors. In our series of works shown by Wang et al. (2020a), we already got the conclusion that the warp is a long lived and not steady structure by adopting all populations and fix the line of node at 5 deg to do the Monte Carlo fitting. It was a relatively simple investigation, but the conclusion was similar to the present work.
In short, to some extent, our results and implications are similar to those obtained by Reylé et al. (2009) Reylé et al. 2009;Momany et al. 2006).  gave the larger negative value of -28 deg. Our results are different from some works corresponding larger values of the warp line of node, but it is well consistent with the value of 5 ± 10 deg used by López-Corredoira et al. (2014) and supporting strongly the works finished by López-Corredoira et al. (2002b). It is intriguing to report here that the red clump stars with 2MASS finished by López-Corredoira et al. (2002b) yielded φ w = −5 ± 5 deg.
Li et al. (2020) also get a different value around 12 deg by using Poggio et al. (2017) simplified model, the difference is caused by that the models and assumptions have much differences, but Li et al. (2020) find that the thick disk population is weaker than the thin disk, which is consistent with our main conclusion.
A three-dimensional map of the Milky Way with the help of classical Cepheid variable stars and the simple model of star formation in the spiral arms was used to reveal the shape of the young stellar disk in . It also mapped the distribution of Cepheid tracers with age for which the range is within a few hundred M yr, much smaller than the range of ages in our sample. Therefore, this is the first time here the warp evolution is followed with complete age sampling of the Milky Way.
In the future, we will go farther than 20−25 kpc of the disk with the state of art of warp model and more accurate and larger sample. We think we can constrain R(kpc) Figure 6. Model and data comparison in six age populations during this work, from the top left one to the bottom right one is corresponding to median value of [1, 3 ,5, 7, 9, 12] Gyr: the black curve is the observed vertical velocity distribution with Galacto-centric radial distance, the blue one is the model profile, the cyan one is the model profile plus 1σ and the green is model profile plus −1σ. The oscillations in the model are due to the fact that the model with sine and cosine function and different values of R have different ranges of φ for most stars.
our galaxy warp better and better, as we mentioned, since there are still relatively large errors in our results. We could also compare stellar warp with gas disk and dust disk warp signals. For the target of this work, we just explore the warp variation with age as an attempt by using a simplified model.

CONCLUSION
In this paper, using LAMOST−Gaia combined red clump giant stars with unprecedented proper motion and age accuracy, we investigate the evolution of warp amplitude, line of nodes, the greatest height and its variation with age. The greatest height of the warp is decreasing with the age and increase with distance in mono-age populations: the younger populations are strongly warped than the old ones. And we also observe the amplitude's temporal evolution and its first derivative with time have a decreasing and increasing pattern respectively. A stable azimuth of line-of-node is −5 degree is shown in this work. Our results are similar to some of recent works, but we use the standard candles with age to quantify warp amplitude evolution.
All these observational results are supporting the warp is not a transient structure, and also implying strongly that warp evolution is non-uniform, long lived, nonsteady structure. We conclude that the warp might be induced by the non-gravitational interaction scenarios: gas infall onto the disc (López-Corredoira et al. 2002b) or magnetic fields (Battaner & Jiménez-Vicente 1998) or similar classes. It might reflect some puzzling evolution of the warp that should be further studied. Both the simple model and data used here can be improved. We need a better warp model and more accurate age measurements to further research this S-like stellar disk. Our work might be of vital importance for us to investigate more properties of the warp and more work will be shown in our series of works.
We would like to thank the anonymous referee for his/her helpful comments. HFW is supported by the built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by National Astronomical Observatories, Chinese Academy of Sciences. This work has also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.