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

ABSTRACT

We select 1052 469 (754 635) thin disc stars from Gaia eDR3 and LAMOST DR7 in the range of Galactocentric radius R (guiding centre radius R_g) from 8 to 11 kpc to investigate the asymmetries between the North and South of the disc mid-plane. More specifically, we analyse the vertical velocity dispersion profiles (⁠|$\sigma _{v_{z}}(z$|⁠)) in different bins of R (R_g) and [Fe/H]. We find troughs in the profiles of |$\sigma _{v_{z}}(z)$| located in both the North (z ∼ 0.7 kpc) and South (z ∼ −0.5 kpc) of the disc 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_g. A similar shift exists in these North–South (NS) asymmetry profiles further divided into different [Fe/H] ranges. The sample binned with R_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.

1 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 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, Gilmore & Gaia-ESO Consortium 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 disc in dynamical equilibrium.

Widrow et al. (2012) found that the disc has a wave-like NS asymmetry in the star counts based on SDSS data, and the disc has a non-zero vertical bulk motion that slowly changes with vertical height. After carefully considering the influence of observational errors and selection biases, Yanny & Gardner (2013) confirmed the existence of an 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 colour of main-sequence stars. Widmark, Widrow & Naik (2022) used Gaia DR3 data to make non-parametric estimates of the number density and the average velocity field of the disc. They found the results significantly deviate from the axisymmetric and NS mirror symmetric models.

Several early studies found asymmetries in the vertical bulk motion of the disc 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 disc is a combination of these two modes (Widrow et al. 2012, 2014; Carlin et al. 2013; Williams et al. 2013; 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 R ∈ 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 z–v_z colour-coded by either the azimuthal velocity (v_ϕ) or the radial velocity (v_R), which indicates that the disc 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 2023).

In modelling our Galaxy, often the vertical velocity dispersion profile (⁠|$\sigma _{v_{z}}(z$|⁠)) 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; Guo et al. 2020; Salomon et al. 2020). Understanding any asymmetries in (⁠|$\sigma _{v_{z}}(z$|⁠)) is essential as this will affect the derived parameters from the modelling (e.g. Banik, Widrow & Dodelson 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 disc mid-plane. 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 characterized 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 disc are considered to be important heating mechanisms (Spitzer & Schwarzschild 1953; Barbanis & Woltjer 1967; Lacey 1984; Aumer, Binney & Schönrich 2016a, b; 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 z–v_z phase space (Antoja et al. 2018; Binney & Schönrich 2018; Tian et al. 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 disc, stars are pulled away from the centre of phase space and produce a high-density region in the disc with the similar phase angle in the z–v_z phase space. For stars oscillating in the disc’s anharmonic vertical gravitational potential, those with different vertical action J_z, namely, different phase space orbits, will have different oscillation frequency Ω_z. This results in a differential rotation in phase space, and the perturbed stars evolve into a phase spiral feature. The Ω_z − J_z 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 z direction. Additionally, the phase spiral winds tighter at smaller radii (Bland-Hawthorn et al. 2019; Wang et al. 2019; Li & Shen 2020; Xu et al. 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, 2023; Tian et al. 2018; Li 2021; Darragh-Ford et al. 2023; Frankel 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 analyse the features in the vertical velocity dispersion profiles (⁠|$\sigma _{v_{z}}-z$|⁠) in different bins of Galactocentric radius R (or guiding centre radius R_g) and metallicity ([Fe/H]). We investigate the variation of the asymmetric vertical velocity dispersion profiles with R (R_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 dependence. 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 (x = −8.34, y = 0, z = 0.027) kpc (Chen et al. 1999; Reid et al. 2014) with a velocity relative to the Local Standard Rest (LSR) as |$(U=9.58, V=10.52, W=7.01)\, \rm 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 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 centre radius.

2 DATA

In order to obtain six dimensional kinematic information, we utilize the precise parallaxes and proper motions in Gaia eDR3 (Gaia Collaboration 2021), and the line-of-sight velocities provided by LAMOST DR7 (Cui et al. 2012; Luo et al. 2015). The cross-matched catalogue between Gaia eDR3 and LAMOST DR7 contains 10 143 169 stars. We exclude the duplicate sources by self-cross-matching and keep the star with highest signal-to-noise ratio in g band (SNR_g).

We apply the following criteria to ensure the data quality:

• parallax_err/parallax < 0.2,

• pm_err/pm < 0.2,

• SNR_g > 20,

where pm is the proper motion defined as |$\sqrt{{\tt pm\_ra}^2+{\tt pm\_dec}^2}$|⁠.

We use [Fe/H] > −0.4 as the condition to select thin disc 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 LSR: |V − V_LSR| < 200 km s⁻¹. 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⁻¹ 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 anticentre direction (Deng et al. 2012; Zhao et al. 2012; Liu, Zhao & Hou 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\rm \, kpc\lt R\lt 11 \rm \, 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 centre radius (R_g), which reflects the orbital properties of stars. We calculate the guiding centre radius using galpy (Bovy 2015) under the Milky Way potential MWPotential2014. The radial criterion we use to define the Sample 2 is |$8\rm \, kpc\lt \mathit{ R}_\mathrm{g}\lt 11\rm \, kpc$|⁠. Fig. 1 shows the R − z and R − ϕ distribution of the full matched LAMOST/Gaia stellar catalogue and of our selected Sample 1 and Sample 2. Our Sample 1 and Sample 2 contain 1052 469 and 754 635 stars, respectively.

The influence of the bias in distance on the phase spiral shape is discussed in Antoja et al. (2023). 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. 2023). 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 Fig. 2, we show the R − z map colour-coded by the median vertical velocity error and median z error. The median of the vertical velocity error is |$\rm \lt 5\, km\, s^{-1}$| and the median of the z error is |$\rm 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-of-sight velocity error is dominant in the vertical direction, and is larger than the error of the precise proper motions.

3 RESULTS

The characteristics of NS asymmetry detected in the vertical velocity dispersion (⁠|$\sigma _{v_{z}}$|⁠) profile give clues to its origin. In this section, we characterize 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.

3.1 Peaks and troughs in |$\sigma _{v_{z}}-z$| profile

We calculate the vertical velocity dispersion |$\sigma _{v_{z}}$| from the standard deviation in v_z. Considering the uncertainty of vertical velocity will lead to an overestimated standard deviation, we correct the median uncertainties of v_z as |$\sigma _{v_{z}}=(\sigma ^2 -\bar{v}_{z,\mathrm{err}}^2)^{1/2}$|⁠. We estimate the uncertainty on |$\sigma _{v_{z}}$| using the bootstrapping method. Black data points in both the left and right panels of Fig. 3 show the vertical velocity dispersion profiles in different ranges of Galactocentric radius R (the left panel) and guiding centre radius R_g (the right panel). Apart from the general trend that the vertical velocity dispersion increases with increasing |z|, the most striking feature is the troughs of the profile near z = 0.6 kpc and z = −0.5 kpc. The features are clearer in the profiles binned with R_g. This is reasonable since the stars have similar orbits in a given R_g bin. The same peak and trough features are found in subsamples of different ranges in radius and metallicity, as shown in Fig. 4. Thus, we rule out the possibility that these peaks and troughs are due to individual stellar associations. Although the shape of the |$\sigma _{v_{z}}$| profile is independent of metallicity, Fig. 4 shows an overall shift in the |$\sigma _{v_{z}}$| 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 z–v_z phase space as a superposition of a smooth background and an asymmetric spiral component. If the spiral intersects with the z-axis at a certain z, there are more stars with v_z 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 z-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 z-axis, which will enhance the trough found in the |$\sigma _{v_{z}}$| profile. We expect that the peak and trough features that we find in the |$\sigma _{v_{z}}-z$| curve thus trace the phase spiral.

3.2 NS asymmetry of |$\sigma _{v_{z}}$|

Since the heights of the intersections between the phase spiral and the z-axis are different in the North and South, the heights of the peaks and troughs in the |$\sigma _{v_{z}}$| profile are different between the North and South, as shown by the blue and red curves in the left (R bins) and right (R_g bins) panels of Fig. 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 R and R_g, and in Section 3.4, bins of [Fe/H]. We find that such comparisons further highlight the characteristics of the NS asymmetry.

Fig. 5 shows the difference between the North and South velocity dispersion profiles, i.e.|$\Delta \sigma _{v_{z}}=\sigma _{v_{z},\mathrm{north}}-\sigma _{v_{z},\mathrm{south}}$| (we refer to this as the NS asymmetry curve of the velocity dispersion profile), binned along R (upper panel) and R_g (lower panel). The NS asymmetry curve is characterized by a wave feature that first shows a peak then trough with increasing |z|. 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 |z|. The wave at smaller R (R_g) curves is shifted to increasingly larger values of |z|. This shift is clearer in the asymmetry curves binned with R_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 mid-plane of the disc is not obvious for different R (R_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 R (R_g) in Fig. 5. The formation of the phase spiral originates from the anharmonicity of the vertical oscillation of stars in the disc (Antoja et al. 2018). The oscillation frequency of stars decreases with the increase of J_z. In other words, in the |$\rm 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 centre is slower. The spiral feature comes from the differential rotation of stars with different J_z. Since the J_z − Ω_z 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 R may have different orbital properties, such that the guiding centre 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 R_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 R_g. Stars with smaller R_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 R_g located in the outer edge, which is well shown in fig. 18 of Li (2021), so that the vertical location of the z-axis intersection is larger. The NS asymmetric curve of stars with smaller R_g will be biased to a relatively larger |z|.

3.3 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 Fig. 6 and the left column of Fig. 7.

We define the density contrast as δρ = n/n_G − 1, where n is number count of stars in the phase space cells in (z, v_z) and n_G is the count smoothed by a Gaussian filter with σ_G = 3 cell width. Fig. 6 shows the phase spirals within three R_g bins colour-coded by δρ. The dashed lines show the locations of the intersection between the spiral and z-axis (same lines plot in the left panels of Fig. 7). The first row of Fig. 7 shows the |$\Delta \sigma _{v_{z}}$| profile in three R_g bins, which is same as the bottom of the Fig. 5. In the second row of Fig. 7, we use δρ_{(N − S)} 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 |$|v_{z}|\lt 10\, \rm km\, s^{-1}$|⁠. The lower panel left column of Fig. 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 |$\sigma _{v_{z}}$| 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 an NS asymmetry (lower panel of Fig. 7). The peaks and valleys of the NS asymmetry wave pattern of |$\sigma _{v_{z}}$| are at similar |z| 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 z − v_z phase space have higher metal abundance. After obtaining a perturbation in velocity, the stars gain energy and enter a phase space orbit with higher J_z. 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 disc with the same J_z. 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 z-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 z − V_z phase space by colour-coding based on metal abundance. Gaia Collaboration (2023) used [M/H] residuals to colour-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 centre radius. In our calculations of the correlations, we exclude the data with |z| > 1.0 kpc due to the large uncertainties. We first build 50 subsamples, each one consisting of a random selection of 60 per cent of the stars from our full sample. We calculate 50 sets of NS asymmetry profiles from different R_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 p-values we calculate using this bootstrapping method does not account for the significant autocorrelation of the signals; this implies that both error bars and p-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 Fig. 7. The correlations between the NS asymmetry of number density and metallicity are low in R_g ∈ (9, 11) kpc, while we find a moderate correlation in R_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 |z|, which can decrease the Pearson correlation coefficient (⁠|$P_\mathrm{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.

3.4 Metallicity dependence of NS asymmetry of |$\sigma _{v_{z}}$| and phase spiral

Fig. 8 shows the NS asymmetry curves of vertical velocity dispersion for Sample 1 (upper row, binned by R) and Sample 2 (lower row, binned by R_g) separated into two [Fe/H] bins. In the top panels of Fig. 8, the wave feature tends to shift to higher |z| with increasing [Fe/H], while the bottom panels slightly show such a tendency in the range of 9 kpc < R_g < 10 kpc (this only slight shift in the R_g-binned samples is due to the stars sharing more similar orbital properties than those binned by R). For stars with different metallicities, their R (R_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], R_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 |z|. Nevertheless, as the samples are constrained to the same R or R_g bins, the shift of peaks between different metallicity bins is not as clear as those of the different radial bins shown in Fig. 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. Fig. 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 |$r=\sqrt{(v_z/55\ \mathrm{km\, s}^{-1})^2+(z/1\ \mathrm{kpc})^2}$|⁠. The phase space is divided into a grid of (θ, r) with an r width of 0.04 kpc. The (θ, r) 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 per cent to calculate the median values and the uncertainties of r. Finally, we convert the uncertainties to the errors in z, v_z.

As shown in the right panel of Fig. 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 R_g (V_ϕ). Stars with poorer metallicity are then related to a shallower vertical potential. This is consistent with the result in Antoja et al. (2023) that the shape of the phase spiral varies with R_g. These results are also consistent with the shift in the wave pattern of the NS asymmetry |$\sigma _{v_{z}}$| curves, because the peaks and troughs of the |$\sigma _{v_{z}}$| profile are actually the intersections between the phase spiral and z-axis.

4 TEST PARTICLE SIMULATIONS

In this section, we use test particles to simulate the evolution of a disturbed thin disc 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 N-body simulations. Although 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 N-body simulation. With these caveats, we proceed to employ the use of test particle simulation.

4.1 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 R in the NS asymmetry of |$\sigma _{v_{z}}$|⁠. The mass and velocity of the intruder are 2 × 10¹⁰ M_⊙ and 300 km s⁻¹, respectively. The impact point is at (x, y) = (20, 0) kpc and the impact angle is 90°. The passage time-scale is T = 66 Myr. The in-plane disturbance velocity δv_|| is simplified as the product of the disturbance time-scale T and the relative acceleration, which is reduced by the acceleration at the Galactic Centre. The vertical disturbance velocity is calculated as δv_⊥ = α δv_||[1 − βsin (θ − θ_Intruder)], 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 disc particles that obey a quasi-isothermal distribution function (Binney 2010; Binney & McMillan 2011). The radial exponential disc 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 kpc < R < 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] that varies with R_g and z to mock a [Fe/H] distribution in z–v_z phase space qualitatively similar to observations (Gaia Collaboration 2023). According to the initial position of the particles in the z–v_z phase space before perturbation, we give a random [Fe/H] value obeying a normal distribution with a mean of μ = −0.1(R_g − 8) − 0.4|z| and a dispersion of σ = 0.1 dex, where R_g and z is in kpc. We integrate the particles under the Galactic potential for 0.6 Gyr until they form a stable [Fe/H]–z distribution, after which we introduce the vertical velocity disturbance.

4.2 NS asymmetries induced by the spiral

Fig. 10 shows the evolution of the test particles in the radial range 9 kpc < R_g < 10 kpc. The top row of Fig. 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 z-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 z 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.

4.3 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 Figs 7 and 11. The right panels of Fig. 7 display the asymmetries for particles with |$\phi \in (-\frac{3\pi }{4},-\frac{\pi }{4})$|⁠, and show small shifts similar to the observations.

The stars entering in a current spatial volume (R − ϕ 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 |$\Delta \sigma _{v_z}$| profiles.

By modelling the multiple impacts of a Sagittarius dwarf-like galaxy through a cold Milky Way-like disc using an N-body simulation, Bland-Hawthorn & Tepper-García (2021) show that the v_ϕ- and v_r-traced phase spiral develops both in the clockwise and anticlockwise direction due to the disc’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 R_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. 2023). The same limitation in ϕ for particles with different R_g will produce different epicyclic phase selection effects. In the range of |$\phi \in (-\frac{3\pi }{4},-\frac{\pi }{4})$|⁠, the spiral is more loosely wound in 8 kpc < R_g < 9 kpc as compared to the more tightly wound spiral in the radial range of 10 kpc < R_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 R_g. This readily implies metallicity dependence of the phase spiral. By using a [Fe/H] label that only depends on R_g, i.e. μ = −0.1(R_g − 8), we reproduce the metallicity dependence of |$\Delta \sigma _{v_z}$| profiles. Fig. 12 shows the metallicity dependence from the simulation. The particles are in the R 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.

5 DISCUSSION

5.1 Impact of selection effects

The primary selection effects in this work are first, those due to the sky coverage of the cross-matched sample, and second, those propagated through from the spectroscopic survey, particularly on stellar colours, magnitudes, and distances. Due to the location, survey strategy, and science goals of LAMOST, there is a strong focus on the anti-Galactic Centre direction and the Northern sky. We inspect the spatial distributions in the R − z (and R_g − z) space and the histograms of R (and R_g) of our samples, normalized by the total stellar numbers in the North and South, respectively. There is only a significant difference in the R distribution of the subsample in the range of 8 kpc < R < 9 kpc due to the survey footprint. The distributions form R_g are quite similar, even in the range of 8 kpc < R_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 v_z direction. Considering a narrow z 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 σ_z profiles. In addition, we also check the density contrast profiles (1D) and maps (2D) in Figs 7 and 9, respectively. Although 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.

5.2 Differences between simulation and observations

In this section, we note the differences in the location of the peaks and troughs of the |$\Delta \sigma _{v_z}$| profiles in the simulation compared to the observations and we consider several explanations. The shift between peaks is larger in the simulation.

The position of the peaks and shifts depend on the orientation of the spirals. The orientation of the spiral involves the initial phase angle of the perturbation, the wrapping speed (or frequency) at a certain radius, and the evolution time. As indicated by recent works (Antoja et al. 2023; Darragh-Ford et al. 2023; Frankel et al. 2023), the evolution time of spirals varies with different angular momenta (R_g). This difference will also change the orientations of the spirals, and thus change the shifts of the peaks.

Self-gravity might change the wrapping frequencies and the orientations of spirals. Darling & Widrow (2019) found that self-gravity can extend the duration of the phase spiral due to the persistent vertical oscillations. After about 1 Gyr of evolution, the phase spiral without self-gravity winds tighter than a spiral with self-gravity. A subsequent theoretical work by Widrow (2023) shows that self-gravity can change the shape of the phase spiral and slow down the phase mixing rate. The evolution time inferred with self-gravity would be shorter than that derived from simulations without self-gravity. The effect that the self-gravity slows down the phase mixing may change the orientation and shape of the phase spirals, and thus reduce the shift in the |$\Delta \sigma _{v_z}$| profiles. However, the extent of such an influence due to self-gravity is still unclear, and may vary with different mechanisms of the phase spiral. The consideration of self-gravity is non-trival and is beyond the current scope of this work.

5.3 Literature comparison

The NS asymmetry in stellar number density was first discovered by Widrow et al. (2012) and later confirmed by more detailed investigations (Yanny & Gardner 2013; Bennett & Bovy 2019). In these works, the number density profiles show excesses at z ∼ −0.4 kpc and z ∼ 0.7 kpc, resulting in the trough and peak of the NS number density difference profile at |z| ∼ 0.4 kpc and |z| ∼ 0.7 kpc, respectively. The locations of the peaks and troughs presented in our results are consistent with these works. An (2019) first found the NS asymmetry in the mean metallicity using the cross-matched sample of SDSS 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 parametrized model to fit the |$\sigma _{v_z}$| 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 |$\sigma _{v_z}$| profiles by this phenomenological spiral model, which also gives an NS density difference profile consistent with the observations. Both studies show the strong connection between the NS asymmetry and the phase spiral.

Investigating the z–v_z 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 disc. They conclude that both the radial metallicity gradient and the age-velocity dispersion relation affect the metallicity dependence. Gaia Collaboration (2023) also found the metallicity dependence and the spiral feature in the vertical phase space colour-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 |$\sigma _{v_z}$| 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 |$\sigma _{v_z}$|⁠, metallicity, and the stellar number density. In addition, we focus more on the shift in the peak locations of |$\Delta \sigma _{v_z}$| 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 z-axis. Guo et al. (2024) directly measured the intersections between the phase snails and z/v_z axes, and utilized an interpolation method to derive the vertical potential for intersections at different radii. The good consistency with previous popular Milky Way potential models demonstrates that the shift in the intersections, i.e. the shift of peak locations in this work, is due to the radial variation of the vertical potential. Therefore, 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.

6 SUMMARY

Using thin disc 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 centre radial range of 8–11 kpc. We find peak and trough features in |$\sigma _{v_z}-z$| profiles that appear at all radial ranges and in subsamples separated by [Fe/H]. We demonstrate that the peaks and troughs in the |$\Delta \sigma _{v_z}$| profiles are indicators of the intersections between the phase spiral and z-axis. The difference in |$\sigma _{v_z}$| between the North and South, i.e. the |$\Delta \sigma _{v_z}-z$| NS asymmetry profile, shows a wave-like pattern. We find strong correlation between the |$\Delta \sigma _{v_z}-z$| 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 |$\Delta \sigma _{v_z}$| profiles binned by Galactocentric radius R and by guiding centre radius R_g. The peak locations of stars with smaller R (R_g) are shifted to larger values of |z|. This trend becomes even more obvious when binning the sample with R_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 metal-rich stars have smaller guiding centre 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 |$\Delta \sigma _{v_z}$| 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 |$\Delta \sigma _{v_z}$| 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 disc 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 R_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. 2023; Darragh-Ford et al. 2023; Frankel 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 (e.g. Guo et al. 2024), as well as the presentation of carefully constructed perturbation modelling.

ACKNOWLEDGEMENTS

We thank Feng Wang, Qiang Yuan, and Zhaoyu Li for helpful discussions and acknowledge the National Natural Science Foundation of China (NSFC) under grant nos. 12373033, 12303023, 11873034, 12011530421, and 12103031, the Department of Science and Technology of Hubei Province for the Outstanding Youth Fund (2019CFA087), the Cultivation Project for LAMOST Scientific Payoff and Research Achievement of CAMS-CAS, and the science research grants from the China Manned Space Project including the CSST Milky Way and Nearby Galaxies Survey on Dust and Extinction Project CMS-CSST-2021-A09 and No. CMS-CSST-2021-A08. RG is supported by Initiative Postdocs Supporting Program (No. BX2021183), funded by China Postdoctoral Science Foundation. This work has 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. Guoshoujing 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.

DATA AVAILABILITY

The data supporting this article will be shared upon reasonable request sent to the corresponding authors.

References