Mapping the Galactic Disk with the LAMOST and Gaia Red Clump Sample. VI. Evidence for the Long-lived Nonsteady Warp of Nongravitational Scenarios

Most spiral galaxies have a warped disk (Sánchez-Saavedra et al. 1990, 2003; Reshetnikov & Combes 1998), although in some spirals, warps cannot be observed due to their low inclination. Like many spiral galaxies in the universe, the warp of the Milky Way (MW) was detected by neutral hydrogen (H i) gas many years ago (Kerr 1957; Bosma 1981; Briggs 1990; Nakanishi & Sofue 2003; Levine et al. 2006b). Dust has also been observed in Marshall et al. (2006). Furthermore, there are several works that measure the 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 three-dimensional (3D) map of the Galactic warp traced by classical Cepheids (Chen et al. 2019; Skowron et al. 2019a), including precession measurements consistent with the Briggs rule (Briggs 1990) and showing us that the amplitude of its 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 Schönrich & Dehnen (2018) and Huang et al. (2018) in the vertical velocity, angular momentum, azimuthal velocity, or guiding center radius parameter space, showing that the Galactic warp is not a transient structure, 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 are 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) and Shen & Sellwood (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 & Binney 1989; Quinn & Binney 1992; Bailin & Steinmetz 2003). Burke (1957) and Weinberg & Blitz (2006) suggested it was caused by the interaction with the Magellanic Clouds, and Bailin (2003) proposed the cause as interaction with Sagittarius (Bailin 2003). Hunter & Toomre (1969) claimed that the mass of the Magellanic Clouds is not enough to explain the observed amplitude of the warp; however, new elements of amplification were introduced by Weinberg & Blitz (2006). Other scenarios, that is to say, perturbations by dwarf satellites (Weinberg 1995; Shen & Sellwood 2006), intergalactic magnetic field 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 & Pfenniger (2004), etc., also appeared. These and other mechanisms were proposed, but some observational evidence may favor some of them, and the different interpretations of the formation of the warp are still being hotly debated. In any case, we know that the kinematical distributions of vertical bulk motions will be contributed by warp asymmetrical structure (Wang et al. 2018a, 2020a); in return, vertical motions can be used to constrain the warp's properties.

Warp in the MW bends upward and downward in the northern and southern hemispheres separately with different amplitude, at least in the gas (Levine et al. 2006a). The amplitude of the warp clearly increases strongly with radius and varies with the azimuth angle (López-Corredoira et al. 2014; Liu et al. 2017a; Li et al. 2020). 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 & Dehnen 2018). 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 an angle with respect to the Galactic plane, with vertical amplitude z_w(R, ϕ) = γ (R − R_⊙)^α$\sin (\phi -\phi$_w), where R is the galactocentric distance and ϕ is the galactocentric azimuth. With the calculation of γ, $\dot{\gamma }$ by assuming ϕ_w and α are constant, their work indicated that the main S-shaped structure of the warp is most likely 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: Brown et al. 2018; Gaia Collaboration: Katz et al. 2018), together with spectroscopically classified red clump stars, e.g., the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST; Cui et al. 2012; Deng et al. 2012; Zhao et al. 2012; Liu et al. 2014). In addition to the unprecedented proper motion, we have the stars' ages 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.'s (2014) simple model and the LAMOST DR4 and Gaia DR2 common red clump stars to investigate the warp amplitude and its first derivative and greatest height variation with different age populations so that we can get a better constraint for γ, $\dot{\gamma }$, and ϕ_w. Hence, we could offer our interpretation of the formation and evolution history of MW warp.

This paper is structured as follows. Section 2 is about how we select our red clump stars, velocity derivation, and vertical velocity distribution in different age populations. The model and method are introduced in Section 3. Our results will be shown in Section 4 and discussions are displayed in Section 5. Finally, we give the conclusions of this work.

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), and Zhao et al. (2012). Now we are entering into phase II. Its fiber is 33, the mean seeing during LAMOST observations is around 3'', and the spatial resolution of LAMOST should be around 5''. The selection functions of LAMOST are almost flat along apparent magnitude; more details can also be found in Carlin et al. (2012), Yuan et al. (2015), and Liu et al. (2017b). We select stars in the LAMOST DR4; its 3461 observed plates and stellar parameters for 6,597,527 spectra are derived in Xiang et al. (2017b). The total sample size is 7,620,612, including stars and galaxies.

The red clump star selection details can be found in Huang et al. (2015b, 2020), and the distance and age are determined by the kernel principal component analysis (KPCA) method, which can be found in more detail in Xiang et al. (2017a, 2017b). According to Huang et al. (2020), the distance uncertainties are around 5%–10%, and the age uncertainties are around 30% and were used quite well in Wang et al. (2019, 2020a, 2020b). Red clump stars are well-known horizontal stars and standard candles, so it is not strange that the distance error of our sample is small. As described in Huang et al. (2020), with the help of positions in the metallicity-dependent effective temperature–surface gravity and color–metallicity diagrams, red clump stars could be selected. Supervised by the high-quality asteroseismology data and KPCA method, ages are determined. Using the properties of intrinsic absolution magnitude, the extinction with star pairs method (Yuan et al. 2013), and recalibration of the Ks absolute magnitudes considering both the metallicity and age dependences, we could acquire the distance with uncertainties of 5%–10%. The Gaia DR2 catalog contains high-precision 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⁻¹ 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 consisting of primary red clump giants (RCGs) with few contaminations of secondary red clumps and red giant branch (RGB) stars, the typical purity and completeness of our primary red clump sample are greater than 80%, and we focus on the kinematics but not the star counts, so the completeness may have a very minor influence on the study. We use the latest sample of this catalog to carefully investigate our scientific target according to the following criteria.

• 1.  
  Samples without parameters such as distance, radial velocity, temperature, and surface gravity are removed.
• 2.  
  Stars located inside $| Z|$ < 1 kpc and 8 kpc < R < 14 kpc are chosen.
• 3.  
  Stars with LAMOST spectroscopic signal-to-noise ratio (S/N) > 20 and age less than 15 Gyr are included.
• 4.  
  [Fe/H] > −1.3 dex.
• 5.  
  V_R = [−150, 150] km s⁻¹, V_θ = [−50, 350] km s⁻¹, and V_Z = [−150, 150] km s⁻¹.

We derive the 3D velocities assuming the location of the Sun is R_⊙ = 8.34 kpc (Reid et al. 2014) and Z_⊙ = 27 pc (Chen et al. 2001). The Tian et al. (2015) solar motion values [U_⊙, V_⊙, W_⊙] = [9.58, 10.52, 7.01] km s⁻¹ and other solar motions (e.g., Huang et al. 2015a) will not change our conclusion at all. The circular speed of the local standard of rest (LSR) is adopted as 238 km s⁻¹ (Schönrich 2012). Cartesian coordinates on the basis of coordinate transformation described in Galpy (Bovy 2015) with LAMOST radial velocity have a precision better than 5 km s⁻¹, which is more convenient and direct than those described one by one in Wang et al. (2018a, 2019, 2020a, 2020b). The vertical angular momentum distributions in per unit mass of the sample associated with the error analysis in the R, Z plane of the Cartesian coordinate system are shown in the left panel of Figure 1. Meanwhile, the projected age and [Fe/H] and their measured error distributions are displayed in the middle and right panels, respectively. It reasonably implies that the angular momentum of stars increases with radial distance in the disk, including the corresponding errors, due to disk total momentum being approximate to vertical momentum. The age distribution for which the stars inside are relatively older and those outside are generally younger supports the inside-out formation of the MW disk. The age error is smaller compared with the age value; though there are probably some possible systematic errors that remain to be ignored, it can still give us a good chance to map the dynamical structures in different age populations.

Zoom In Zoom Out Reset image size

As shown in the bottom right panel of Figure 1, with Galactic radial distance increasing, the metallicity has a negative trend, and the error is also becoming larger with a 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) and recent Gaia 3D kinematics works (Gaia Collaboration: Katz et al. 2018; López-Corredoira & Sylos Labini 2019). In this work, we focus on the range of radial distance R = [8 14] kpc, and vertical height Z = [−1 1] kpc; the bins with a minimum number of every pixel containing five stars are shown in Figure 1.

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 Figure 2. As shown in each panel, the vertical velocity increases with radial distance in different mono-age populations from 0 to 10 km s⁻¹ at 1 Gyr. For the others, it is from 0 to 6−8 km s⁻¹, but for the oldest, it is even less; the maximum is around 5 km s⁻¹. They are definitely reflecting the warp signals. And the vertical velocities clearly increase for most age bins around 6−8 km s⁻¹ 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⁻¹. Although there are some oscillations due to the Poissonian noise, as we will mention again in Section 4, in order to get more points to show 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 the 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. During this work, the size and number of the bins are not constant for the different age sample; errors are very different in different populations. We have calculated the number to see their variation labeled in Figure 2. There are many more young stars than old stars, which would explain that the Poisson noise seems smaller at ages less than 5 Gyr. For the population at 9 Gyr, the number is very small, so the large error bar could be caused by Poisson noise. The drop at large R around 12.5 kpc for the last age bin might be caused by the Sun is not being on the line-of-nodes thus cause some stars not to move toward anticenter with different directions possibly so we think it is significant.

Zoom In Zoom Out Reset image size

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, we can see many more details in López-Corredoira et al. (2014). The modeling process is displayed as follows:

$$\begin{eqnarray}&&{V}_{Z}={\rm{\Omega }}(R,z^{\prime} =z-{z}_{w})\sin [{i}_{w}(R)]\cos (\phi -{\phi }_{w})+\dot{{z}_{w}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{z}_{w}(R,\phi )=R\tan [{i}_{w}(R)]\sin (\phi -{\phi }_{w}),\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&{z}_{w}(R\gt {R}_{\odot },\phi ={\phi }_{w}+\pi /2)\approx \gamma {\left(R-{R}_{\odot }\right)}^{\alpha },\end{eqnarray} \tag{ 3 }$$

and a variable line of nodes that has no extremely slow precession to do fitting (i.e., $\dot{{\phi }_{w}}\ll \dot{\gamma }$) and changed with the shape of the warp is 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 of these formulae and assumptions, we derive that, for low angles i_w(R), the low height disk warp model can be simplified reasonably as

$$\begin{eqnarray}\begin{array}{rcl}{V}_{Z}(R\gt {R}_{\odot }) & \approx & \displaystyle \frac{{\left(R-{R}_{\odot }\right)}^{\alpha }}{R}[\gamma {{\rm{\Omega }}}_{\mathrm{LSR}}\cos (\phi -{\phi }_{w})\\ & & +\,\dot{\gamma }R\sin (\phi -{\phi }_{w})],\end{array}\end{eqnarray} \tag{ 4 }$$

where ϕ_w is the azimuth of the line of nodes, γ is the amplitude of the warp, and $\dot{\gamma }$ describes the warp amplitude evolution, that is to say, −d(γ)/d(age). We assume the exponent α = 2 (kpc⁻¹) (López-Corredoira et al. 2014) and set the line of nodes as ϕ_w (deg) as a free parameter, ϕ_w = variable degree (in the literature the values are between −28°, −5°, +15°, and 18°; Chen et al. 2019; López-Corredoira et al. 2002b; Momany et al. 2006; Reylé et al. 2009; Skowron et al. 2019a). The constant α = 1 (no units; Reylé et al. 2009) is also tested in our studies, and we checked that the conclusions derived from the fitting results are stable and not much affected by the value of α. Here we just use the data with R ≥ 8.34 kpc beyond the Sun to get the best-fitting value based on a Markov Chain Monte Carlo (MCMC) simulation provided by EMCEE (Foreman-Mackey et al. 2013). The model and method are used maturely by our series of works in Wang et al. (2020a) by fitting all populations. With the help of the carefully selected sample, the present model, and the popular MCMC method, we could obtain the likelihood distribution of the vertical velocity profile for fitting as

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{l}{ \mathcal L }(\{{V}_{{\rm{obs}}}({R}_{i}| Z,{\rm{Age}})\}| \gamma ,\dot{\gamma },{\phi }_{w})=\displaystyle \prod _{i}\\ {\rm{\exp }}\left[-\displaystyle \frac{1}{2}\left({V}_{{\rm{obs}}}({R}_{i}| Z,{\rm{Age}}\right)-{\left.{V}_{{\rm{model}}}({R}_{i}| Z,{\rm{Age}},\gamma ,\dot{\gamma },{\phi }_{w}\right)}^{2}\right].\end{array}\end{array}\end{eqnarray*} \tag{ 5 }$$

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 in Figure 2 to constrain the warp for this work. Please notice that we set the parameter as $\dot{\gamma }$ = −d(γ)/d(age) by adopting a joint likelihood with 12 parameters to do simultaneous fitting of all age bins. In order to get the convergent parameters and save computer time, we just choose a relatively smaller sampling size in our simulation; the MCMC size is 50*12*1000, and the step is 500. 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.

4.1. 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. 2018a, 2020a), which implies that clearly vertical upward motions can be used to reveal the properties of the warp, such as warp amplitude, precession rate, and so on. The fitting results of the work in different age bins, by fitting all of the age bins simultaneously with α = 2 in the model, are displayed in Figure 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.

Zoom In Zoom Out Reset image size

4.2. Warp Parameter γ, $\dot{\gamma }$, and ϕ_w Evolution with Age

In Figure 4, the amplitude γ (${\mathrm{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; the 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 recalculated again by bootstrap process. We can see that there is a variation of the γ (kpc⁻¹) with age with a relatively large 1 − σ error; all of these values are decreasing with age. Correspondingly, the bottom left panel is the warp amplitude derivative variation $\dot{\gamma }$ (kpc⁻¹ Gyr⁻¹) distribution with age; it has a variable increasing trend. These values are different from zero, implying that the warp always exists but is not a stationary structure ($\dot{\gamma }$ (kpc⁻¹ Gyr⁻¹) ∼ 0), and there is a clear difference of populations existing. Moreover, there is also a stable feature for the azimuth of the line of nodes in all populations; the value is almost fixed at about −5° for the distribution displayed in the top right panel. The variation is very small due to the relative larger error and is shown as a flat line, and some theoretical studies also support that the line of nodes is expected to be straight within R $\leqslant$ 4.5 disk scale lengths (Shen & Sellwood 2006; Bland-Hawthorn & Gerhard 2016). The variation of ϕ_w due to precession is too small to be detected in all of the warp models that are also mentioned in López-Corredoira et al. (2002a), Dubinski & Chakrabarty (2009), Jeon et al. (2009), and Poggio et al. (2020). Please note that there are only minor oscillations in the figure when we zoom in, and there might be some intriguing physics in it.

Zoom In Zoom Out Reset image size

We are considering the variation of ϕ_w with time negligible. We cannot distinguish in the vertical motion maps which is the dominant factor in V_Z, $\dot{\gamma }$ or $\dot{{\phi }_{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 & Chakrabarty 2009; Jeon et al. 2009). Moreover, the fits depend on the model for warp shape. In Poggio et al. (2020), they are also taking into account the variation of the amplitude of warp. They have a long-lived precessing model with a large value by assuming a single value and no radial motions with the help of an average of four geometrical simplified models. We suggest that the calculation of the precession by Poggio et al. (2020) with the very young population of Chen et al. (2019) is assuming it is similar for old and young population, also implied in Extended Data Figure 3 for similar precession values of four models, and the model of warp they have used with Gaia DR2 is not good because it does not reach high R, a good description of fits results of the old population in comparison with the young population warp is displayed in Chrobáková et al. (2020), which is also showing younger populations like Cepheids are worse than whole populations. There is a possibility that when they get too-high values of the vertical motions without precession, they have to introduce a too-high precession to compensate for it. We are skeptical about the validity of the results of Poggio et al. (2020), and they need to be investigated more.

As a natural product, the $\dot{\gamma }$/γ is shown in the bottom right panel, which is different from zero and increases with age. The figures derived from the MCMC simulations for all of these age bins could give us a reliable estimation for the parameters of γ (kpc⁻¹) and ϕ_w (deg) thanks to the peak of the maximum likelihood. Note that γ has units of kpc⁻¹ for our adopted value of α = 2; it has no units when we set α = 1.

4.3. Z_w: Greatest Height of Warp for Younger Ages

In Figure 5, we show the distribution of the greatest height of the 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 height along with distance in all age bins. We also see that the younger populations are higher than the old ones. Note that there are six age bins, but two lines are overlapped. For the bottom panel, we can clearly see that there is clear decreasing pattern for all median heights in all age populations. It is consistent with the results in Figure 4, meaning that the warp appears to be a long-lived nonstationary structure, and, more importantly, there is a clear difference for different populations again.

Zoom In Zoom Out Reset image size

4.4. The Comparison of Model and Sample

In order to check the robustness of the derived results, we have finished the comparison of the model of Equation (4) and data in six age populations during this work, as shown in Figure 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 solid line with error bars is the observed vertical velocity distribution with galactocentric radial distance in each panel; the blue line is the model profile with the Monte Carlo fitting, the cyan line is the model profile plus 1σ, and the green line is the 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. The few mismatches shown in Figure 6 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 toward anticenter with different directions; extinctions in some regions are not perfect.

Zoom In Zoom Out Reset image size

In order to try to test the Poissonian noise to mislead us for the conclusion, we just compare our model with data consisting of at least 100 stars in each bin, which has small differences for the pattern shown in Figures 2 and 6. The overall trend is better except that fewer boundary points deviate, which cannot change our conclusion at all. So we suggest that all features mentioned in the previous paragraphs are robust and intriguing. It seems that there is a systematical peak in the model at R = 12 kpc of Figure 6; we actually have a test for it, suggesting that our model with sine and cosine functions and different values of R has different ranges of ϕ for most stars. If we have a simple Monte Carlo simulation, there are some clear oscillations and peaks located around 12 kpc, so it is expected and possibly caused by our model distribution properties.

Therefore, in short summary, we think that these stable features are real, and this observational evidence strongly supports that the warp is a long-lived and nonsteady S-shaped structure; it also implies that the warp evolution is relatively uniform and there is a clear difference for different age bins. So far, we only have a few points spanning 8–14 kpc with relatively large and nonperfect age errors; it would be worthwhile for us to further investigate this structure in more detail with the help of a sample consisting of a larger distance range, more stars, and a more accurate age. Some scenarios are given in the next section.

Vertical nonaxisymmetries and wavelike density patterns are found in the solar neighborhood and the outer disk (Widrow et al. 2012; Carlin et al. 2013; Williams et al. 2013; Xu et al. 2015; Pearl et al. 2017; Carrillo et al. 2018, 2019; Wang et al. 2018a, 2018b, 2018c). As mentioned and implied in Wang et al. (2020a), many mechanisms including warp might be coupled together under the same complex distribution function to cause the complicated galactoseismology signals. Scenarios for producing these structures, such as warp dynamics, minor mergers, or interactions with nearby dwarfs or satellites (Gómez et al. 2013; D'Onghia et al. 2016; Laporte et al. 2018, 2019), and the effects of even lower-mass dark matter subhalos 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 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 anticenter region are discussed in Liu et al. (2017a), where it was found that the vertical bulk motion of younger red clump stars is significantly larger than that of the older ones, which is consistent with our results. We have a nonsteady 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 favor 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 a larger warp amplitude.

According to Skowron et al. (2019b), there are mainly two classes of warp formation mechanisms. One is the warp formed by gravitational interactions, for example, with satellite galaxies or a misaligned dark matter halo. The other one is nongravitational mechanisms, e.g., accretion of intergalactic gas onto the disk (López-Corredoira et al. 2002a) or interactions with intergalactic magnetic fields. Nongravitational mechanisms such as magnetic fields or hydrodynamical pressure from infalling gas would act on the gas and only affect the young stars (Guijarro et al. 2010; Sellwood 2013); thus, we should see that all signals of younger ones are stronger than old ones, and we do get that clearly. Therefore, we think the gravitational scenario should not be the reason, at least for this tracer.

The young population traces the gas, whereas the 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 the gas and not the stars. Our current results support that the warp might be contributed by the nongravitational 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 on both the position and kinematics of the Galactic warp has been observed by Amôres et al. (2017) and Romero-Gómez et al. (2019) and confirmed for our results in Figures 2 and 5. We could clearly see that there are some vertical velocity differences and the greatest height difference in different age bins. Amôres et al. (2017) thought that the warp shapes and amplitude of the northern part are stronger than the southern part, and the north–south asymmetry was also presented in Reylé et al. (2009) and Momany et al. (2006), which we cannot test here. The reason is that LAMOST mainly covers the Galactic anticenter in the north, and we do not have many stars in the southern sky; it is impossible to discuss this difference, but we plan to work it out in the future. Romero-Gómez et al. (2019) showed that the amplitude of OB stars corresponding to younger populations is weaker than the RGB stars corresponding to older ones by calculating the onset radius and height of the warp, thus suggesting that the warping disk of the older population is more pronounced or stronger. This 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) used 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, Skowron et al. (2019a) showed that their results, by using Cepheids similar to OB stars, are in contradiction to Romero-Gómez et al. (2019) in that the Cepheids' warp height is similar to RGB stars.

We also notice that the recent results of Chen et al. (2019) for 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 nonsteady structure by adopting all populations and fixing the line of nodes at 5° 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), Amôres et al. (2017), Chen et al. (2019), Liu et al. (2017a), Chrobáková et al. (2020) and others, although we have some differences with Romero-Gómez et al. (2019), Poggio et al. (2018), and Skowron et al. (2019a), with some possible reasons. Further work needs to be done to clarify these differences.

In Chen et al. (2019), their line of nodes is around 18° with the help of lines of nodes in different radial bins. In other literature, the values are between −5° and +15° (López-Corredoira et al. 2002b; Momany et al. 2006; Reylé et al. 2009; Li et al. 2020). Skowron et al. (2019a) gave the larger negative value of −28°. Our results are different from some works with corresponding larger values of the warp line of nodes, but they are well consistent with the value of 5° ± 10° used by López-Corredoira et al. (2014) and strongly support 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 ${\phi }_{w}=-5^\circ \pm 5^\circ$.

Li et al. (2020) also got a different value around 12° by using the Poggio et al. (2017) simplified model; the difference is because the models and assumptions have many differences, but Li et al. (2020) found that the warp signal of the thick disk population is weaker than the thin disk, which is consistent with our main conclusion.

A 3D map of the MW 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 Skowron et al. (2019a). It also mapped the distribution of Cepheid tracers with age for which the range is within a few hundred Myr, much smaller than the range of ages in our sample. Therefore, this is the first time the warp evolution is followed with complete age sampling of the MW.

In the future, we will go farther than 20−25 kpc of the disk with the state of the art of the warp model and a more accurate and larger sample. We think we can constrain our galaxy warp better and better, as we mentioned, since there are still relatively large errors in our results. We could also compare the stellar warp with gas disk and dust disk warp signals. Moreover, we could use [Fe/H] and [α/Fe] as population indicators to see more evolution features of the warp. For the target of this work, we just explore the warp variation with age as an attempt to use a simplified model.

In this paper, using LAMOST and Gaia combined red clump giant stars with unprecedented proper motion and age accuracy, we investigate the evolution of the warp amplitude, line of nodes, and greatest height and its variation with age. The greatest height of the warp is decreasing with age and increasing with distance in mono-age populations: the younger populations are more strongly warped than the old ones. And we also observe that the amplitude's temporal evolution and its first derivative with time have a decreasing and increasing pattern, respectively. A stable azimuth of the line of nodes is −5°, as shown in this work. Our results are similar to those of some recent works, but we use the standard candles with age to quantify warp amplitude evolution.

All of these observational results support that the warp is not a transient structure and also imply strongly that warp evolution is a nonuniform, long-lived, nonsteady structure. We conclude that the warp might be induced by the nongravitational interaction scenarios: gas infall onto the disk (López-Corredoira et al. 2002b), 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 thank the anonymous referee for the helpful comments. H.F.W. is supported by the LAMOST Fellow project, National Key Basic R&D Program of China via 2019YFA0405500 and funded by China Postdoctoral Science Foundation via grant 2019M653504, Yunnan province postdoctoral Directed Culture Foundation and the Cultivation Project for LAMOST Scientific Payoff and Research Achievement of CAMS-CAS. Y.H. acknowledges the National Natural Science Foundation of China U1531244,11833006, 11811530289, U1731108, 11803029, and 11903027 and the Yunnan University grant No. C176220100006 and C176220100007. M.L.C. was supported by grant PGC-2018-102249-B-100 of the Spanish Ministry of Economy and Competitiveness. J.L.C. acknowledges support from the U.S. National Science Foundation via grant AST-1816196. H.W.Z. is supported by the National Natural Science Foundation of China under grant No. 11973001 and supported by National Key R&D Program of China No. 2019YFA0405500. H.F.W. is fighting for the plan "Mapping the Milky Way Disk Population Structures and Galactoseismology (MWDPSG) with large sky surveys" in order to establish a theoretical framework in the future to unify the global picture of the disk structures and origins with a possible comprehensive distribution function. We pay our respects to elders, colleagues, and others for comments and suggestions; thanks to all of them. The Guo Shou Jing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope, LAMOST) is a National Major Scientific Project 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 the 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.