North-South asymmetries in the Galactic thin disk associated with the vertical phase spiral as seen using LAMOST-Gaia stars

We select 1,052,469 (754,635) thin disk stars from {\it Gaia} eDR3 and LAMOST DR7 in the range of Galactocentric radius $R$ (guiding center radius $R_\mathrm{g}$) from 8 to 11\,kpc to investigate the asymmetries between the North and South of the disk midplane. More specifically we analyze the vertical velocity dispersion profiles ($\sigma_{v_{z}}(z$)) in different bins of $R$ ($R_\mathrm{g}$) and $[\mathrm{Fe/H}]$. We find troughs in the profiles of $\sigma_{v_{z}}(z)$ located in both the North ($z \sim 0.7$\,kpc) and South ($z \sim -0.5$\,kpc) of the disk at all radial and chemical bins studied. The difference between the Northern and Southern vertical velocity dispersion profiles ($\Delta\sigma_{v_{z}}(|z|)$) shows a shift between curves of different $R$ and $R_\mathrm{g}$. A similar shift exists in these NS asymmetry profiles further divided into different $[\mathrm{Fe/H}]$ ranges. The sample binned with $R_\mathrm{g}$ more clearly displays the features in the velocity dispersion profiles. The shift in the peaks of the $\Delta\sigma_{v_{z}}$ profiles and the variation in the phase spiral shape binned by metallicity indicate the variation of the vertical potential profiles and the radial metallicity gradient. The wave-like signal in NS asymmetry of $\sigma_{v_{z}}(z)$ largely originates from phase spiral; while the NS asymmetry profiles of [Fe/H] only display a weak wave-like feature near solar radius. We perform a test particle simulation to qualitatively reproduce the observed results. A quantitative explanation of the NS asymmetry in the metallicity profile needs careful consideration of the spiral shape and the perturbation model, and we leave this for future work.


INTRODUCTION
The vertical density and kinematics of the Milky Way play an important role in understanding the structure and dynamics of the Milky Way as a whole.Data from large-scale sky surveys, e.g., Large Sky Area Multi-Object Fiber Spectroscopic Telescope survey (LAMOST, Cui et al. 2012;Yan et al. 2022), Gaia (Gaia Collaboration et al. 2021), Sloan Digital Sky Survey (SDSS, Almeida et al. 2023), RAdial Velocity Experiment (RAVE, Steinmetz et al. 2006), Gaia-ESO (Gilmore et al. 2012;Randich et al. 2013), GALactic Archaeology with HERMES (GALAH, De Silva et al. 2015;Buder et al. 2021), are vital to reveal the complex details of the vertical structure of the Milky Way; these data break the simple assumption of a North-South (NS) symmetric Milky Way disk in dynamical equilibrium.Widrow et al. (2012) found that the disk has a wave-like NS asymmetry in the star counts based on SDSS data, and the disk has a non-zero vertical bulk motion that slowly changes with vertical ★ Email: hjtian@lamost.org, guorui20@sjtu.edu.cn,sarahbird@ctgu.edu.cn,liuchao@nao.cas.cn.height.After carefully considering the influence of observational errors and selection biases, Yanny & Gardner (2013) confirmed the existence of a NS asymmetry in the stellar number counts.With the large sample size and high-precision parallaxes from Gaia DR2, Bennett & Bovy (2019) found that the asymmetry in number density is independent of the color of main sequence stars.Widmark et al. (2022) used Gaia DR3 data to make non-parametric estimates of the number density and the average velocity field of the disk.They found the results significantly deviate from the axisymmetric and NS mirror symmetric models.
Several early studies finding asymmetries in the vertical bulk motion of the disk include.Many studies of the vertical bulk motion have found that there are two modes, i.e, the breathing mode and bending mode.The vertical bulk motion of the disk is a combination of these two modes (Widrow et al. 2012;Carlin et al. 2013;Williams et al. 2013;Widrow et al. 2014;Wang et al. 2018;López-Corredoira et al. 2020).Ding et al. (2021) selected K giant stars from LAMOST and studied the vertical kinematic structure in the Galactocentric radial range of  ∈ 5 − 15 kpc and up to 3 kpc vertically from the Galactic plane.They found that the vertical velocity dispersion in the South is larger than that in the North in the radial range of 7 − 13 kpc.Antoja et al. (2018) found a spiral feature in the phase space of  −   color-coded by either the azimuthal velocity (  ) or the radial velocity (  ), which indicates that the disk is undergoing a phase mixing process possibly initiated by a vertical perturbation.
Besides the asymmetry in kinematics and density, the vertical distribution of metallicity in the Solar Neighbourhood is found to show a wave-like asymmetry, which coincides with the previously known asymmetry in the stellar number density distribution (An 2019).In addition, the phase spiral structures in the density and kinematic distributions were found to be linked to the metallicity ([M/H]) variations (Gaia Collaboration et al. 2022).
In modelling our Galaxy, often the vertical velocity dispersion profile (   ()) is needed, such as in using the vertical Jeans equation to estimate the vertical force and local surface density and dark matter density (e.g., Widrow et al. 2012;Hagen & Helmi 2018;Salomon et al. 2020;Guo et al. 2020).Understanding any asymmetries in (   ()) is essential as this will affect the derived parameters from the modelling (e.g., Banik et al. 2017;Sivertsson et al. 2018;Haines et al. 2019).Salomon et al. (2020)'s Jeans analysis shows that velocity dispersion is more perturbed in the North than in the South of the disk midplane.They additionally find a connection between the NS asymmetry in number density and the vertical phase spiral.Guo et al. (2020) found the velocity dispersion curve shows a plateau feature at different heights in the North and South.These plateaux result in different dark matter density estimates between the North and South.
The vertical velocity dispersion distribution can reflect the history of dynamical heating, as characterised by the age-velocity dispersion relationship.The small scale irregular gravitational potential generated by giant molecular clouds and the transient spiral arm structures in the disk are considered to be important heating mechanisms (Spitzer & Schwarzschild 1953;Lacey 1984;Barbanis & Woltjer 1967;Aumer et al. 2016a;Aumer et al. 2016b;Ting & Rix 2019).External perturbation is also an important source of heating (e.g., Gómez et al. 2012;Grand et al. 2016).Such perturbation can cause the bending modes of vertical bulk motion, and initiate the formation of a phase spiral in the  −   phase space (Antoja et al. 2018;Tian et al. 2018;Binney & Schönrich 2018;Bland-Hawthorn et al. 2019;Laporte et al. 2019;Bland-Hawthorn & Tepper-García 2021).
The formation mechanism and the properties of the vertical phase spiral have been studied in great depth.Antoja et al. (2018) and Binney & Schönrich (2018) elaborated upon the likely origin of the phase spiral.As a high-speed dwarf galaxy impacts the disk, stars are pulled away from the center of phase space and produce a high-density region in the disk with the similar phase angle in the  −   phase space.For stars oscillating in the disk's anharmonic vertical gravitational potential, those with different vertical action   , namely, different phase space orbits, will have different oscillation frequency Ω  .This results in a differential rotation in phase space, and the perturbed stars evolve into a phase spiral feature.The Ω  −   distribution depends on the vertical gravitational potential, which means that stars with similar orbital properties will produce a similar phase spiral shape.The gravitational potential will become shallower with increasing radius, resulting in different spiral shapes at different radii, and more elongated spirals in the  direction.Additionally, the phase spiral winds tighter at smaller radii (Bland-Hawthorn et al. 2019;Wang et al. 2019;Xu et al. 2020;Li & Shen 2020).Under a given gravitational potential, the phase spiral shape will evolve with time.Thus, one can infer the impact time from the shape of the phase spiral (Antoja et al. 2018;Tian et al. 2018;Li 2021;Antoja et al. 2022;Frankel et al. 2022;Darragh-Ford et al. 2023).Tian et al. (2018) constrained the vertical perturbation to starting no later than 0.5 Gyr ago.The aforementioned observed characteristics of the phase spiral are in good agreement with theory.Guo et al. (2022) investigated the possible connection between the NS asymmetry of the velocity dispersion and the vertical phase spiral found by Antoja et al. (2018).Guo et al. (2022) established a model through the superposition of an equilibrium background and a phase spiral component to quantitatively explain the asymmetries in the number density and velocity dispersion profiles.
In this work, we compile a large sample cross-matched from the recent LAMOST and Gaia data releases, to specifically analyze the features in the vertical velocity dispersion profiles (   − ) in different bins of Galactocentric radius  (or guiding center radius  g ) and metallicity ([Fe/H]).We investigate the variation of the asymmetric vertical velocity dispersion profiles with  ( g ) and [Fe/H], and finally utilize a test particle simulation to confirm the connection between the NS asymmetry of the velocity dispersion and the vertical phase spiral found in the observations.This paper is arranged as follows.In Section 2, we introduce the data and sample selection criteria.In Section 3, we present the results from our data analyses of the NS asymmetries, including the NS asymmetry of the vertical velocity dispersion and its metallicity dependency.We present how these results are theoretically related to the phase spiral and provide supporting observational evidence.We reproduce these results qualitatively through test particle simulation as presented in Section 4, and detail the necessary factors in order to reproduce the characteristics similar to the observations.Finally, we provide a discussion and summarize in Sections 5 and 6, respectively.
Throughout the paper, we adopt the position of the Sun as ( = −8.34, = 0,  = 0.027) kpc (Chen et al. 1999;Reid et al. 2014) with a velocity relative to the Local Standard Rest (LSR) as ( = 9.58,  = 10.52, = 7.01) km s −1 (Tian et al. 2015).These adopted values are well consistent with those widely recognized works (e.g., Schönrich 2012;Bland-Hawthorn & Gerhard 2016;GRAVITY Collaboration et al. 2019).In principle, the solar motion with respect to the LSR does not affect the velocity dispersion, the position of the Sun just slightly affects the calculation of the guiding center radius.

DATA
In order to obtain six dimensional kinematic information, we utilize the precise parallaxes and proper motions in Gaia eDR3 (Gaia Collaboration et al. 2021), and the line-of-sight velocities provided by LAMOST DR7 (Cui et al. 2012;Luo et al. 2015).The cross-matched catalog between Gaia eDR3 and LAMOST DR7 contains 10,143,169 stars.We exclude the duplicate sources by self-crossmatching and keep the star with highest signal-to-noise ratio in  band (SNR  ).
We apply the following criteria to ensure the data quality: where pm is the proper motion defined as √︁ pm_ra 2 + pm_dec 2 .We use [Fe/H] > −0.4 as the condition to select thin disk stars, where [Fe/H] is that provided by LAMOST DR7.Note that our sample includes both main sequence and giant stars.In order to exclude halo stars, we adopt a kinematic cut in the absolute velocity relative to the local standard of rest (LRS): | − LSR | < 200 km s −1 .The stars remaining after this step are defined as the preliminary sample (named as Sample 0).Note that the line-of-sight velocity  from LAMOST is systematically smaller than that from Gaia.We add 5.7 km s −1 to our LAMOST line-of-sight velocities as a bias correction (Tian et al. 2015).
A key focus of LAMOST is to observe in the Galactic anticenter direction (Deng et al. 2012;Zhao et al. 2012;Liu et al. 2015), resulting in a dearth of stars at small Galactocentric radii.To ensure a sufficient number of stars for binning the data, we avoid the inner Galactic regions.We select stars located in the radial range of 8 <  < 11 kpc and the azimuthal angle range of −10 • <  < 10 • in the Galactocentric cylindrical coordinates as our first sample (named as Sample 1 throughout the paper).The second sample (named as Sample 2) is defined similarly and only differs in the definition of the radius: we select stars using the guiding center radius ( g ), which reflects the orbital properties of stars.We calculate the guiding center radius using galpy (Bovy 2015) under the Milky Way potential MWPotential2014.The radial criterion we use to define the Sample 2 is 8 <  g < 11 kpc.Figure 1 shows the  −  and  −  distribution of the full matched LAMOST/Gaia stellar catalog and of our selected Sample 1 and Sample 2. Our Sample 1 and Sample 2 contain 1,052,469 and 754,635 stars, respectively.
The influence of the bias in distance on the phase spiral shape is discussed in Antoja et al. (2022).Distances derived from the inverse of parallaxes suffer from an asymmetric systematic bias, which is then transferred to the calculation of coordinates and velocities.Failure to account for this effect will cause the phase spiral to shrink (see the appendix of Antoja et al. 2022).To avoid this effect, we therefore adopt the distances from Bailer-Jones et al. ( 2021).We convert the distances, line-of-sight velocities, and proper motions to cylindrical coordinates and corresponding uncertainties by adopting the method of Johnson & Soderblom (1987).In Figure 2 we show the  −  map color-coded by the median vertical velocity error and median  error.The median of the vertical velocity error is < 5 km s −1 and the median of the  error is 0.06 kpc in the majority of regions studied, which attests to the good quality of our data set.We note that the error increases with Galactic latitudes.This is because the line-ofsight velocity error is dominant in the vertical direction, and is larger than the error of the precise proper motions.

RESULTS
The characteristics of NS asymmetry detected in the vertical velocity dispersion (   ) profile give clues to its origin.In this section we characterise the vertical velocity dispersion profile and look for clues to a common origin with the vertical phase spiral.We elaborate upon the new features that we find within the NS asymmetry profile.Incorporating the origin of the phase spiral, we provide a qualitative explanation for the NS vertical velocity dispersion asymmetry and its detailed characteristics.Finally, we show the connection between NS asymmetries (including stellar density and metallicity) with the phase spiral.

Peaks and troughs in 𝜎 𝑣 𝑧 − 𝑧 profile
We calculate the vertical velocity dispersion    from the standard deviation in   .Considering the uncertainty of vertical velocity will lead to an overestimated standard deviation, we correct the median uncertainties of   as    = ( 2 − v2 ,err ) 1/2 .We estimate the uncertainty on    using the bootstrapping method.Black data points in both the left and right panels of Figure 3 show the vertical velocity dispersion profiles in different ranges of Galactocentric radius  (the left panel) and guiding center radius  g (the right panel).Apart from the general trend that the vertical velocity dispersion increases with increasing ||, the most striking feature is the troughs of the profile near  = 0.6 kpc and  = −0.5 kpc.The features are clearer in the profiles binned with  g .This is reasonable since the stars have similar orbits in a given  g bin.The same peak and trough features are The data are binned using a linear gridding with width of 0.05 kpc.Right panel: Same as the left panel, except using Sample 2, thus the data are binned using the guiding center radius  g .Compared to the results in the left panels, the features (e.g., the troughs near  = 0.6 kpc and  = −0.5 kpc) are more striking in the profiles binned with  g shown in the right panels .found in subsamples of different ranges in radius and metallicity, as shown in Figure 4. Thus, we rule out the possibility that these peaks and troughs are due to individual stellar associations.Although the shape of the    profile is independent of metallicity, Figure 4 shows an overall shift in the    profile with metallicity.
The phenomenological model proposed by Guo et al. (2022) is able to explain the trough features of the vertical velocity dispersion curve by associating them with the components of the asymmetric phase spiral.Guo et al. (2022) simply regards the distribution of stars in the  −   phase space as a superposition of a smooth background and an asymmetric spiral component.If the spiral intersects with the  axis at a certain , there are more stars with   close to zero velocity, which will reduce the velocity dispersion.For vertical regions away from the intersections, the velocity dispersion will increase.Therefore, the phase space spiral will cause a peak in the vertical velocity dispersion profile at the intersection between the phase spiral and the  axis.In addition, stars in the process of vertical oscillation are more often at high vertical height and low velocity states, thus the phase space spiral component will be enriched near the intersection of the phase space orbit and the  axis, which will enhance the trough found in the    profile.We expect that the peak and trough features that we find in the    −  curve thus trace the phase spiral.

NS asymmetry of 𝜎 𝑣 𝑧
Since the heights of the intersections between the phase spiral and the  axis are different in the North and South, the heights of the peaks and troughs in the    profile are different between the North and South, as shown by the blue and red curves in the left ( bins) and right ( g bins) panels of Figure 3.This results in the NS asymmetry of vertical velocity dispersion.We now investigate the difference between the North and South velocity dispersion profiles by comparing our samples in bins of different  and  g , and in Section 3.4, bins of [Fe/H].We find that such comparisons further highlight the characteristics of the NS asymmetry.
Figure 5 shows the difference between the North and South velocity dispersion profiles, i.e., Δ   =    ,north −    ,south (we refer to this as the NS asymmetry curve of the velocity dispersion profile), binned along  (upper panel) and  g (lower panel).The NS asymmetry curve is characterized by a wave feature which first shows a peak then trough with increasing ||.We find that the same wave feature of the NS asymmetry curve appears in all three radial bins investigated, although the feature is slightly shifted along ||.The wave at smaller  ( g ) curves is shifted to increasingly larger values of ||.This shift is clearer in the asymmetry curves binned with  g .The intensity of the NS asymmetry produced by the phase spiral component is affected by the ratio of the stellar number in the spiral components and the background.The NS asymmetry of velocity dispersion near the midplane of the disk is not obvious for different  ( g ) bins probably because the ratio of the spiral component to background is small.
If we assume that the asymmetry of the vertical velocity dispersion is caused by the phase spiral, and combine this with the formation mechanism and characteristics of the phase spiral, we are able to qualitatively explain the shifts seen in the curves at different  ( g ) in Figure 5.The formation of the phase spiral originates from the  anharmonicity of the vertical oscillation of stars in the disk (Antoja et al. 2018).The oscillation frequency of stars decreases with the increase of   .In other words, in the z − v z phase space, the angular velocity of a star near the origin is faster, while the angular velocity of a star far from the center is slower.The spiral feature comes from the differential rotation of stars with different   .Since the   − Ω  depends on the vertical gravitational potential profile, the evolution of the phase space spiral fundamentally depends on the vertical gravitational potential experienced by these stars.Therefore the phase spiral shapes at different Galactic radii are different due to the different vertical potential (Wang et al. 2018).However, stars at the same  may have different orbital properties, such that the guiding center radius is a better way to distinguish phase spirals with different intrinsic shapes (Li & Shen 2020;Hunt et al. 2021).Separating the samples according to  g to get phase spiral components with the same intrinsic shape reduces the morphological ambiguity brought about by mixing.Thus, the shift in the NS asymmetry curves are clearer when samples are binned by  g .Stars with smaller  g experience steeper vertical gravitational potential, so their oscillation frequency is generally higher and they tend to rotate with larger phase angles at the same time interval.Therefore, the spiral consists of stars with small  g located in the outer edge, which is well shown in Figure 18 of Li (2021), so that the vertical location of the -axis intersection is larger.The NS asymmetric curve of stars with smaller  g will be biased to a relatively larger ||.

NS asymmetries and their connection with the phase spiral
In the previous sections, we have detailed the NS asymmetry of the vertical velocity dispersion observed in our sample and have made inferences under the assumption of a direct connection to the vertical phase space spiral.Here we explore the NS asymmetries in number density and [Fe/H] and their connection with the phase spiral.Our results are summarized in Figure 6 and the left column of Figure 7.
We define the density contrast as  = /  − 1, where  is number count of stars in the phase space cells in (,   ) and   is the count smoothed by a Gaussian filter with   = 3 cell width.Figure 6 shows the phase spirals within three  g bins color-coded by .The dashed lines show the locations of the intersection between the spiral and  axis (same lines plot in the left panels of Figure 7).The first row of Figure 7 shows the Δ   profile in three  g bins, which is same as the bottom of the Figure 5.In the second row of Figure 7, we use  (  −) to represent the NS asymmetry in number density, which is the difference of  between the North and South for a stripe in phase space with |  | < 10 km s −1 .The lower panel left column of Figure 7 shows the NS asymmetry curve for [Fe/H].The peaks in the NS asymmetry of the vertical velocity dispersion roughly correspond to the trough in the NS asymmetry of number density and metallicity.As indicated by the dashed lines, the common shift shown in the NS asymmetry profiles of both    and  is strong evidence of a common origin due to the phase spiral.
Similar to the vertical velocity dispersion distribution and density distribution, the vertical distribution of metallicity abundance also shows a NS asymmetry (lower panel of Figure 7).The peaks and valleys of the NS asymmetry wave pattern of    are at similar || oppositely fluctuating for  and [Fe/H].The oppositely fluctuating asymmetry in metallicity were first reported by An (2019).At the similar location where the wave feature in the vertical velocity dispersion NS asymmetry curve shows a peak, there is a valley in the NS difference of both density and [Fe/H].This phenomenon is a natural consequence of the model describing the asymmetry as a product of the vertical phase spirals excited by disturbances.The phase spiral component is composed of stars that have been disturbed by the gravitational potential of a perturber.Stars in the central part of the  −   phase space have higher metal abundance.After obtaining a perturbation in velocity, the stars gain energy and enter a phase space orbit with higher   .As the orbits of these excited stars evolve, these stars become the main components of the phase spiral and they have higher [Fe/H] than the undisturbed background components of the disk with the same   .Here we ignore the effect caused by the radial disturbance and only consider the change of the metal abundance distribution caused by the vertical velocity excitation.As mentioned earlier in Section 3.1, vertically oscillating stars stay longer at low velocity states, thus the number of spiral components is higher near the  axis in vertical phase space.As a result, the increase of median metal abundance is more obvious.
The spiral also can be displayed in the  −   phase space by color-coding based on metal abundance.Gaia Collaboration et al. ( 2022) used [M/H] residuals to color-code the vertical phase space distribution.They also suggest that this NS asymmetry feature is possibly connected with the phase spiral.
If the wave signal in the NS asymmetry of number density originates from the phase spiral as we assume in this study, the correlation between the NS asymmetry of number density and other asymmetries will reflect how much the phase spiral is able to account for this NS asymmetry signal.To inspect this correlation, we calculate the Pearson correlation between the NS asymmetry of density and velocity dispersion, and the Pearson correlation between the NS asymmetry of number density and metallicity.As this statistic measures the linear correlation coefficient, such correlation does not precisely represent the relation between the NS asymmetry profiles, but is enough to allow us to make a qualitative judgement.
Table 1 shows the correlation between the NS asymmetry curves in different ranges of guiding center radius.In our calculations of the correlations, we exclude the data with || > 1.0 kpc due to the large uncertainties.We firstly build 50 subsamples, each one consisting of a random selection of 60% of the stars from our full sample.We calculate 50 sets of NS asymmetry profiles from different  g bins using the same process as in previous sections.Then we calculate the Pearson correlation between the NS asymmetry of number density and metallicity from each set of subsamples.Finally we adopt the 50th percentile as our result and 16th and 84th percentiles as our error bar.The error bars and -values we calculate using this bootstrapping method does not account for the significant autocorrelation of the signals; this implies that both error bars and -values are underestimated.
The high correlation between the NS asymmetry of number density and vertical velocity dispersion shows their strong connection.The strong connection is also supported by the close agreement of their peak/trough locations, as seen by the dashed lines of Figure 7.The correlations between the NS asymmetry of number density and metallicity are low in  g ∈ (9, 11) kpc, while we find a moderate correlation in  g ∈ (8, 9) kpc.On the one hand, the wake-like signal may be smoothed out by the uncertainties on [Fe/H].This could weaken the peaks and troughs at large radial range.On the other hand, there is a trend that the NS difference varies with ||, which can decrease the Pearson correlation coefficient ( corr_coeff ).Our analysis with Pearson's correlation statistic supports our assumption that the wave-like signal is caused by the phase spiral, but is only detected near the solar radius.

Metallicity dependence of NS asymmetry of 𝜎 𝑣 𝑧 and phase spiral
Figure 8 shows the NS asymmetry curves of vertical velocity dispersion for Sample 1 (upper row, binned by ) and Sample 2 (lower row, binned by  g ) separated into two [Fe/H] bins.In the top panels of Figure 8, the wave feature tends to shift to higher || with increasing [Fe/H], while the bottom panels slightly show such a tendency in the range of 9 <  g < 10 kpc (this only slight shift in the  g -binned samples is due to the stars sharing more similar orbital properties than those binned by ).For stars with different metallicities, their  ( g ) are generally different.The orbital properties of the stars affect the shape of the spiral, and thus lead to the differences in the NS asymmetric curves.For stars with higher [Fe/H],  g is generally smaller due to the negative radial metallicity gradient.The vertical potential they experience is steeper, resulting in a phase spiral located at the outer region of phase space.Therefore, the wave pattern in the NS asymmetry curves of the velocity dispersion shifts to larger ||.Nevertheless, as the samples are constrained to the same  or   bins, the shift of peaks between different metallicity bins is not as clear as those of the different radial bins shown in Figure 5.
Using the preliminary sample (i.e., Sample 0 defined in Section 2) with an additional criterion in azimuth (−10 • <  < 10 • ), we investigate how the shape of the phase spiral changes with metallicity.Figure 9 shows the spirals within two metallicity bins.To obtain the shape of the phase spiral, we evenly divide each spiral into 20 sectors, and depict the shape of phase spiral by the density peaks of each of the 20 sectors.We normalize the polar coordinates of the phase space by  = √︃ (  /55 km s −1 ) 2 + (/1 kpc) 2 .The phase space is divided into a grid of (, ) with an  width of 0.04 kpc.The (, ) with the maximum mean  in each sector is adopted as the coordinates of the density peak.We randomly resample 50 times with a ratio of 70% to calculate the median values and the uncertainties of .Finally, we convert the uncertainties to the errors in ,   .
As shown in the right panel of Figure 9, the spiral consisting of metal-rich stars (blue curve) is generally located on the outward side of the metal-poor spiral (green curve).The metal-rich stars tend to dominate the outer edges along the phase spiral curve.This is a reflection of the radial metallicity gradient, as the metallicity bin can be regarded as a bin in  g (  ).Stars with poorer metallicity are then related to a shallower vertical potential.This is consistent with the result in Antoja et al. (2022) that the shape of the phase spiral varies with  g .These results are also consistent with the shift in the wave pattern of the NS asymmetry    curves, because the peaks and troughs of the    profile are actually the intersections between the phase spiral and  axis.

TEST PARTICLE SIMULATIONS
In this section, we use test particles to simulate the evolution of a disturbed thin disk in the gravitational potential of the Milky Way, we aim to verify whether the phase spiral generated by external disturbances could produce the observed NS asymmetry in velocity dispersion and [Fe/H].We note that simple test particle simulations cannot reproduce the self gravity effect such as in -body simulations.Though the results are not quantitative fits to the observations, they are sufficient enough to produce qualitatively comparable results.For our initial simple simulations, using test particles is much more time and cost efficient compared to a full -body simulation.With these caveats, we proceed to employ the use of test particle simulation.

Models adopted in the simulations
We mimic the impact of an intruder by adopting the perturbation model under the impulse approximation given by Binney & Schönrich (2018).Our goal of using this simulation is to reproduce the observed variation with  in the NS asymmetry of    .The mass and velocity of the intruder are 2 × 10 10 M ⊙ and 300 km s −1 , respectively.The impact point is at (, ) = (20, 0) kpc and the impact angle is 90 • .The passage time scale is  = 66 Myr.The in-plane disturbance velocity  | | is simplified as the product of the disturbance time scale  and the relative acceleration, which is reduced by the acceleration at the Galactic centre.The vertical disturbance velocity is calculated as where we set  = 1 and  = 1.
In our simulations, we adopt model I of Irrgang et al. (2013) for the gravitational potential of the Milky Way and use AGAMA to generate thin disk particles that obey a quasi-isothermal distribution function (Binney 2010;Binney & McMillan 2011).The radial exponential disk scale length is set as 3.7 kpc, and the vertical scale height is set to 0.3 kpc (Binney & Piffl 2015;Bland-Hawthorn & Gerhard 2016).We generate approximately 1.5 million test particles in the range of 5 <  < 15 kpc and implement orbital integration with galpy.
In order to mock the influence of the velocity perturbation on the distribution of metallicity, the test particles are labelled with a random [Fe/H] which varies with  g and  to mock a [Fe/H] distribution in  −   phase space qualitatively similar to observations (Gaia Collaboration et al. 2022).According to the initial position of the particles in the  −   phase space before perturbation, we give a random [Fe/H] value obeying a normal distribution with a mean of  = −0.1(g − 8) − 0.4|| and a dispersion of  = 0.1 dex, where  g and  is in kpc.We integrate the particles under the Galactic potential for 0.6 Gyr until they form a stable [Fe/H] −  distribution, after which we introduce the vertical velocity disturbance.

NS asymmetries induced by the spiral
Figure 10 shows the evolution of the test particles in the radial range 9 <  g < 10 kpc.The top row of Figure 10 shows the contrast density map in phase space.The corresponding vertical velocity dispersion distribution and median [Fe/H] distribution are shown in the middle and bottom rows, respectively.We find peak and trough features similar to the observations.The intersection of the spiral and the  axis corresponds to the decreases (troughs) in the velocity dispersion as indicated by the grey vertical dashed lines.The vertical velocity dispersion increases significantly at the  region where the vertical velocity of the phase spiral component is high.We find a feature of the median metallicity profiles: the median [Fe/H] shows a bump feature at the height corresponding to the trough in the vertical velocity dispersion.Over time, the phase spiral becomes blurred, and the trough in the velocity dispersion curve and bump in the median metallicity profile become weaker.Note that we do not add background stars, thus the intensity of the NS asymmetries are enhanced in our simulations.Nevertheless, the qualitative result is not influenced.

Shifts of the NS asymmetries found in the simulation
We inspect the phase spirals from the simulation at the evolution time of 0.5 Gyr, with the results shown in Figure 7 and 11.The right panels of Figure 7 display the asymmetries for particles with  ∈ (− 3  4 , −  4 ), and show small shifts similar to the observations.The stars entering in a current spatial volume ( −  volume) come from different locations of the Galactic plane, where they have in the past received quite different vertical velocity kicks at the time of the perturbation.The effect of mixing will make the phase spiral more disperse and thus change the peak locations and amplitudes of the Δ   profiles.
By modelling the multiple impacts of a Sagittarius dwarf-like galaxy through a cold Milky Way-like disk using an -body simulation, Bland-Hawthorn & Tepper-García (2021) show that the   -and   -traced phase spiral develops both in the clockwise and anticlockwise direction due to the disk's shear.Our test particle simulation induces a velocity kick and, at an evolution time of 0.5 Gyr, shows that the density-traced phase spiral emerges only in specific regions of azimuth.By binning along  g , our simulation shows that the shapes of the phase spirals wind tighter when the azimuth varies anticlockwise.Such varying shape of the phase spiral with azimuth has been shown in the Milky Way using the Gaia DR3 data (Antoja et al. 2022).The same limitation in  for particles with different  g will produce different epicyclic phase selection effects.In the range of  ∈ (− 3  4 , −  4 ), the spiral is more loosely wound in 8 <  g < 9 kpc as compared to the more tightly wound spiral in the radial range of 10 <  g < 11 kpc.The small shifts in the observations may partly result from these effects.Binney & Schönrich (2018) first point out that the evolution of the phase spiral depends on different vertical frequencies that are different along  g .This readily implies metallicity dependence of the phase spiral.By using a [Fe/H] label that only depends on  g , i.e.  = −0.1(g − 8), we reproduce the metallicity dependence of Δ   profiles.Figure 12 shows the metallicity dependence from the simulation.The particles are in the  range of 8−9 kpc at an evolution time of 0.5 Gyr.The result from the simulation shows a small shift similar to the observations, mainly due to the mixture effect from a radially and azimuthally dependent vertical perturbation.

Impact of selection effects
The primary selection effects in this work are firstly, those due to the sky coverage of the cross-matched sample, and secondly, those propagated through from the spectroscopic survey, particularly on stellar colors, magnitudes and distances.Due to the location, survey strategy, and science goals of LAMOST, there is a strong focus on the anti-Galactic-Center direction and the Northern sky.We inspect the spatial distributions in the  −  (and  g − ) space and the histograms of  (and  g ) of our samples, normalized by the total stellar numbers in the North and South respectively.There is only a significant difference in the  distribution of the subsample in the range of 8 <  < 9 kpc due to the survey footprint.The distributions form g are quite similar, even in the range of 8 <  g < 9 kpc.From this we conclude that such selection effects do not affect the qualitative conclusions in our work.
The complicated selection effects due to the spectral observations will significantly influence the calculation of the real number density profiles.However, in this work, we focus on the velocity dispersion profiles to study the variation in the shape of the phase spiral.The selection effects are trivial in the   direction.Considering a narrow  bin, we conclude the selection effects in this small bin are similar, and also are similar between the North and South.Thus, the selection effects do not greatly influence our results which are derived from the   profiles.In addition, we also check the density contrast profiles (1D) and maps (2D) in Figures 7 and 9, respectively.Though selection effects are clearly present, the Gaussian kernel smoothing helps alleviate such effects to some extent.The connection between the NS asymmetry and the variation in the shapes of the phase spirals are still clear in these two figures.Thus, we conclude that the selection effects do not influence our results and qualitative conclusions.and Gaia DR2.The difference in the mean metallicity displays a wave-like oscillation, which coincides with the NS asymmetry in the stellar number density.An (2019) proposed that the NS asymmetries are induced by the phase spiral.Later, Guo et al. ( 2022) used a parameterized model to fit the    profiles of G/K dwarfs in the Solar Neighbourhood.The model is a superposition of a phase spiral similar to the observations and a smooth background assumed to be in dynamical equilibrium.They quantitatively explained the NS asymmetry in the    profiles by this phenomenological spiral model, which also gives a NS density difference profile consistent with the observations.Both studies show the strong connection between the NS asymmetry and the phase spiral.
Investigating the  −   phase space, Bland-Hawthorn et al. ( 2019) dissected the sample according to the  abundance and the metallicity.They found that the phase spiral shows a clear trend in metallicity: the inner spiral is stronger for metal-rich stars, while the outer spiral is stronger for metal-poor stars, particularly in the -poor disk.They conclude that both the radial metallicity gradient and the age-velocity dispersion relation affect the metallicity dependence.Gaia Collaboration et al. ( 2022) also found the metallicity dependence and the spiral feature in the vertical phase space color-coded by the residual metallicity, which they found to be especially obvious in the more distant radial bins.
Following the work of Guo et al. (2022), we specifically investigate the NS asymmetry in    and the metallicity dependence utilizing the LAMOST data cross-matched with Gaia eDR3.We confirm again the connection between the NS asymmetry and the phase spiral through the coincidence in the NS asymmetry of    , metallicity and the stellar number density.In addition, we focus more on the shift in the peak locations of Δ   profiles.We relate this with the variation of the vertical potential at different radii, as the peak locations are actually the intersections between the phase snail and  axis.We consider the metallicity dependence of the shift of the peak locations and the variation in the shapes of the phase spirals as due to the radial metallicity gradient and the variation of the vertical potential.The peaks and trough features provide an additional way to quantitatively measure the spiral intersections and thus the spiral shape.This will be helpful in relating the spiral shape with a certain vertical potential profile.

SUMMARY
Using thin disk stars selected from LAMOST and Gaia for which we have full 6D information, as well as [Fe/H], we inspect the NS asymmetry of the vertical velocity dispersion in the Galactocentric radial range and guiding center radial range of 8 − 11 kpc.We find peak and trough features in    −  profiles that appear at all radial ranges and in subsamples separated by [Fe/H].We demonstrate that the peaks and troughs in the Δ   profiles are indicators of the intersections between the phase spiral and  axis.The difference in    between the North and South, i.e., the Δ   −  NS asymmetry profile, shows a wave-like pattern.We find strong correlation between the Δ   −  NS asymmetry and number density NS asymmetry, which supports their origin from the phase spiral.On the other hand, we find a weak wave-like signal that is related to the phase spiral from the [Fe/H] NS asymmetry profile near the solar radius.
We also find a variation of the peak locations in Δ   profiles binned by Galactocentric radius  and by guiding center radius  g .The peak locations of stars with smaller  ( g ) are shifted to larger values of ||.This trend becomes even more obvious when binning the sample with  g .This shift in the peak locations reflects the variation in the shape of the phase spiral, which intrinsically indicates the vertical potential at different radii.This shift also exists in the subsamples separated by metallicity, reflecting the influence of the radial metallicity gradient.We separate the total sample into three subsamples according only to metallicity, and roughly measure the shapes of the three phase spirals.The spiral of the metal-rich sample is located on the outward side along the metal-poor spiral.The metalrich stars have smaller guiding center radii and thus a steeper vertical potential, resulting in the spiral lying along the outward side of the metal-poor phase space spiral.The result from the 2D spiral shape is consistent with the results from the 1D Δ   profiles.We perform a test particle simulation similar to that of Binney & Schönrich (2018) in order to qualitatively recover our results found in the observations.Comparable to the observations, we find the wave-like oscillations in the Δ   profiles and the shift in the peak locations for the simulated samples in different radial and chemical ranges.The vertical velocity perturbation can excite the metal-rich stars to higher regions from the disk and produce a metallicity bump relative to the background, which then results in the NS asymmetry in metallicity.The radial and azimuthal variations of the velocity kick and azimuth limitation produce a different epicyclic selection effect for the samples separated into different bins of  g .This selection effect along with the mixing of the stars are likely the main causes of the small amplitude and shift seen in the observations.Our results from the simulation are only qualitatively similar to the observations.Differences remain between the simulation and observations.The spiral shape is influenced by the initial phase angle after the perturbation, the wrapping speed at each radius, and the evolution time since the perturbation.These are found to vary with the angular momentum and the orbital properties of the stars as shown in recent works (Antoja et al. 2022;Frankel et al. 2022;Darragh-Ford et al. 2023).A quantitative explanation of the NS asymmetry in metallicity needs more complicated modelling and quantitative consideration of the relation between the spiral shape and the vertical potential.This will be analysed in future work (Guo et al. in preparation), as well as the presentation of carefully constructed perturbation modelling.

Figure 1 .
Figure 1.Distributions of our samples in Galactocentric cylindrical coordinates.The first and second panels on the left show the distributions of Sample 1 and Sample 2 in  − , respectively.The third and fourth panels show the distributions of Sample 1 and Sample 2 in  − , respectively.The grey markers represent the catalog only after implementing cuts for quality.The red and blue markers represent Sample 1 and Sample 2, respectively.The data selection is detailed in Section 2.

Figure 3 .
Figure 3. Left panel: Comparison of vertical velocity dispersion   profiles in different radial  ranges (as indicated in the top-left corner of in each row) for Sample 1.The left column of the panel shows the vertical velocity dispersion versus .Each data point is estimated by a fixed number  of stars sorted by ;  is indicated in the bottom left of each row.The right column of the panel shows the vertical velocity dispersion profiles separately for both North (blue) and South (red) of the disk midplane binned by |  |.The data are binned using a linear gridding with width of 0.05 kpc.Right panel: Same as the left panel, except using Sample 2, thus the data are binned using the guiding center radius  g .Compared to the results in the left panels, the features (e.g., the troughs near  = 0.6 kpc and  = −0.5 kpc) are more striking in the profiles binned with  g shown in the right panels

Figure 4 .
Figure 4. Vertical velocity dispersion profiles for three different guiding center radial bins of Sample 2. The metallicity range is indicated in the legend of each panel, along with the sample size  per bin.

Figure 5 .
Figure 5. NS   asymmetry profiles for Sample 1 (in different Galactocentric radial ranges, upper panel) and Sample 2 (in different guiding center radial ranges, lower panel), where Δ  =   ,north −   ,south .The green, blue, and black lines represent stars in different radial ranges as indicated in top left corner of each panel.

Figure 6 .
Figure 6.Phase spirals for stars with different  g ranges color-coded by  for our Sample 2. The dashed lines indicate the intersections between the phase spirals and  axis.

Figure 7 .
Figure 7. NS difference profiles from our Sample 2 (left column) compared to our test particle simulation (right column).The solid lines show profiles with different  g ranges, as indicated in the legend.The vertical dashed lines indicate the intersections between the phase spirals with different  g and the  axis, which are the same for Sample 2 as indicated in Figure 6.First row: NS difference profiles of vertical velocity dispersion, i.e.Δ  .Second row: NS difference profiles of number density ( ), for stars within |  | < 10 km s −1 .Third row: NS difference profiles of median [Fe/H].The bins along |  | are 0.05 kpc.For the realistic simulation on the right, we allow the sample to evolve for 0.5 Gyr, with an azimuthal range of [ −3 4 , −1 4  ].

Figure 10 .
Figure10.Results from the test particle simulation.Top row: number density distribution in phase space.Middle row: vertical velocity dispersion profile with error bars.Bottom row: median [Fe/H] profile.The grey lines represent the stable [Fe/H] profiles before we add the vertical velocity kick.The blue lines correspond to the median [Fe/H] profiles after implementing the velocity kick.The first, second, and third columns correspond to the evolution times at 0.4 Gyr, 0.45 Gyr, and 0.5 Gyr, respectively.The vertical dashed lines highlight the coincidence between the phase spiral, the peaks and troughs in the vertical velocity dispersion profiles and the bumps in the [Fe/H] profiles.The radial range of the sample is 9 <  g < 10 kpc and the azimuthal range is  ∈ (− 3  4 , −  4 ).

Figure 12 .
Figure 12.Dependency of the   NS asymmetry with [Fe/H] from simulation after a running time of 0.5 Gyr.The particles from the Galactocentric radial range of 8 ≤  ≤ 9 kpc.We adopt an azimuthal criterion as indicated in the top of the panel.The metallicity [Fe/H] probed by each curve is indicated in the legend.

Table 1 .
Pearson's correlation between NS asymmetry curves in different radial ranges.The uncertainties are estimated by the technique of bootstrapping.8 <  g /kpc < 9 9 <  g /kpc < 10 10 <  g /kpc < 11Note: The error bars and -values we calculate using a bootstrapping method that does not account for the significant autocorrelation of the signals, implying that both error bars and -values are underestimated.