Vertical Phase Mixing across the Galactic Disk

By combining the {\it LAMOST} and {\it Gaia} data, we investigate the vertical phase mixing across the Galactic disk. Our results confirm the existence of the phase space snail shells (or phase spirals) from 6 to 12 kpc. We find that grouping stars by the guiding radius ($R_{g}$), instead of the present radius ($R$) further enhances the snail shell signal in the following aspects: (1) clarity of the snail shell shape is increased; (2) more wraps of the snail shell can be seen; (3) the phase space is less affected by the observational bias, i.e., the lack of stars close to the disk mid-plane due to extinction; (4) the phase space snail shell is amplified in greater radial ranges. The quantitatively measured snail shell shapes are similar, except that the $R_{g}$-based snail shells show more wraps with better contrast. These lines of evidence lead to the conclusion that the guiding radius (angular momentum) is a fundamental parameter tracing the phase space snail shell across the Galactic disk. Results of our test particle simulations with impulse approximation verify that particles grouped according to the guiding radius reveal more prominent and well-defined, sharper snail shell features in the number density phase space. As a guideline for future vertical phase mixing study, it is recommended to use the guiding radius with additional constraints on the orbital hotness (ellipticity) to improve the clarity of the phase space snail shell in the Galactic disk.


INTRODUCTION
The Milky Way, as a massive pure disk galaxy, is not in dynamical equilibrium. The current status of the Galactic disk is mainly shaped via both the internal and external mechanisms together. Resonances from the bar and spiral arms have been shown to induce observed fine structures in the velocity phase space (Dehnen 2000;Fux 2001;Antoja et al. 2009Antoja et al. , 2011Quillen et al. 2011;Hunt & Bovy 2018) and the large scale bulk motions in the Galactic disk (Siebert et al. 2011(Siebert et al. , 2012Carlin et al. 2013;Debattista 2014;Faure et al. 2014;Sun et al. 2015;Monari et al. 2015Monari et al. , 2016Tian et al. 2017;Wang et al. 2018a,b). The external perturbations from satellite galaxies or sub-halos can generate warps, flares, vertical density asymmetries or high order velocity modes in the Galactic disk, such as the bending and breathing vertical motions (Hunter & Toomre 1969;Quinn et al. 1993;Kazantzidis et al. 2008;Purcell et al. 2011;Gómez et al. 2013;D'Onghia et al. 2016;Laporte et al. 2018a,b).
External perturbations also excite vertical phase mixing of stars in the Galactic disk. As first shown in Antoja et al. (2018) with Gaia DR2, stars near the solar radius exhibit a snail shell feature in the Z − V Z phase space, representing the ongoing vertical phase mixing that probably happened 300 − 900 Myr ago. The snail shell forms as a result of anharmonic oscillation of stars in the vertical direction; the vertical oscillation frequency (Ω Z ) is anti-correlated with the vertical action (J Z ) and the angular momentum (L Z ), i.e., V φ for stars in the solar neighborhood (Binney & Schönrich 2018). The phase space snail shell can be recognized for stars in a wide range of ages (Tian et al. 2018;Laporte et al. 2019), and dif-  (Bland-Hawthorn et al. 2019). By dissecting stars in the V R −V φ phase space into distinct arches, Li & Shen (2020) found clear snail shells only in the stars on dynamically colder orbits, i.e., stars with |V φ − V LSR | ≤ 30 km/s (or J R < 0.04 kpc 2 /Myr). The hotter orbits, on the other hand, may have phase wrapped away already due to the much larger radial excursions to facilitate faster phase mixing. The absence of the clear snail shells on hotter orbits suggests that the vertical perturbation occurred at least 500 Myr ago. However, different opinions also exist for the origin of the vertical phase mixing process (see Khoperskov et al. 2019;Bennett & Bovy 2020).
The phase space snail shell has also been found beyond the solar neighborhood in the Milky Way disk. Based on the Gaia DR2 data around the solar neighborhood (within 1 kpc), slightly different phase space snail shells are found at different radius, suggesting the snail shell as a signature of corrugated waves propagating through the disk (Bland-Hawthorn et al. 2019). Laporte et al. (2019) later explored larger radial ranges from 6.6 to 10 kpc, and confirmed the existence of the snail shell at these locations. Similar results can also be seen in Wang et al. (2019) based on the LAMOST and Gaia DR2 data from 6.34 to 12.34 kpc. These studies arrive at similar results; the snail shell at larger radius is more elongated (squashed) along the Z (V Z ) direction and less wound, due to the weaker vertical restoring force and longer orbital time period at larger radius. Recently, Xu et al. (2020) further link the Z −V Z phase space distributions at different radius with the azimuthal bulk motions across the Galactic disk out to 15 kpc, indicating tight connection between the external perturbation and the disk kinematics.
For a volume limited sample, the angular momentum (L Z ), i.e., the guiding radius (R g ), may help to better reveal the global kinematic patterns. For example, based on a local sample near the solar neighborhood, the Galactic warp was identified with the positive correlation with the wave-like  Fig. 2). For the main sample (LAMOST DR6 ∪ Gaia RVS), after crossing-match with the Bayesian distance catalog from Bailer-Jones et al. (2018) and removing stars with large velocity uncertainties, the final sample contains 11,350,423 stars. pattern between the mean vertical velocity and the guiding radius (Schönrich & Dehnen 2018;Huang et al. 2018). Recently, Friske & Schönrich (2019) found that the dependence of the mean radial motion on the angular momentum (guiding radius) also shows wave-like behavior. The shape of the V R − R g profile is similar to, but much stronger than the V R − R profile; patterns in V R − R profile have probably been washed out by excursions of stars around their guiding centers (Friske & Schönrich 2019). Khanna et al. (2019) found that the snail shell shape remains the same at different radius (within 1 kpc from the Sun) for stars with the same L Z , but varies at different L Z . They suggested that the angular momentum could be a more robust quantity in characterizing the snail shell.
Instead of mapping the stellar kinematics and tracing the phase space maps in the traditional spatial volume, grouping stars according to the guiding radius may be another promising way to visualize the morphological or kinematic structures in the disk (e.g., Khoperskov et al. 2020) (but see Hunt et al. 2020 for a different point of view). Aiming to probe the phase space snail shell in a large radial range in the Galactic disk with good number statistics, we combine the Gaia radial velocity sample (RVS) and the LAMOST DR6 sample together with all the stars in both samples included, i.e., LAMOST DR6 ∪ Gaia RVS. Although the LAMOST survey does not sample the sky as uniformly as the Gaia survey, it has good spatial coverage in the Galactic Anti-Center direction with radial velocity information for fainter stars at larger distances, compensating the relatively brighter Gaia RVS sample. With the combined large sample, we can investigate the spatial variation of the snail shell with angular momentum in detail. 3 2. SAMPLE LAMOST DR6 provides reliable radial velocity values for 5,843,107 stars, focusing on the Galactic Anti-Center direction (Zhao et al. 2012;Deng et al. 2012;Liu et al. 2014). The Gaia RVS sample contains 7,225,631 stars with accurate position and velocity information (Gaia Collaboration et al. 2018b). To achieve good statistics, we utilize all the stars to form a main sample of 12,256,045 stars (LAMOST DR6 ∪ Gaia RVS). Fig. 1 summarizes the number of stars in the two samples.  Gaia RVS samples (bottom) in the X − Y plane (left column), the R − Z plane (middle column), and the φ − Z plane (right column). Compared to the relatively brighter Gaia RVS sample, LAMOST sample has its advantage of more faint stars at relatively larger distances.
There are 750,134 stars in common between LAMOST DR6 and Gaia RVS samples (LAMOST DR6 ∩ Gaia RVS). Fig. 2 compares the line-of-sight velocities observed in Gaia (V Gaia LOS ) and LAMOST (V LAMOST LOS ) for these common stars. As shown in the left panel, the agreement between the velocities are quite good. However, there exists a small systematic offset, with the V LAMOST LOS slightly smaller than V Gaia LOS . The distribution of the velocity difference ( is shown in the right panel. The peak position of the distribution is marked with the vertical dashed line. The line-of-sight velocity observed in LAMOST is systematically lower than the Gaia radial velocity by 4.75 km/s. 4 In the following analysis, for these common stars, only the Gaia measured velocities are adopted. For the rest of the stars in the LAMOST DR6 sample, we have added 4.75 km/s to the line-of-sight velocity to compensate this effect. We adopt the Bayesian distance from Bailer-Jones et al. 0.027) kpc as the Sun position (Reid et al. 2014). The local standard of rest (LSR) circular velocity V LSR is set to 240 km/s (Reid et al. 2014). Here the peculiar velocities of the Sun with respect to LSR are set to (U ⊙ ,V ⊙ ,W ⊙ ) = (11.1, 12.24, 7.25) km/s (Schönrich 2012). 5 The spatial distributions of the LAMOST and Gaia RVS samples are shown in Fig. 3. The Gaia RVS sample distributions are quite symmetric in the Y , Z and φ directions. 6 On the other hand, the LAMOST survey mainly covers the Galactic Anti-Center direction and φ > 0 in the positive Y axis. Since the limiting magnitude is deeper in the LAMOST spectroscopic survey, it also extends further in the vertical direction than the Gaia survey.
We further remove stars with relative distance uncertainties larger than 25%, velocity uncertainties larger than 50 km/s, and absolute velocities (|V R |, |V φ −V LSR |, and |V Z |) larger than 400 km/s. The final sample contains 11,350,423 stars in total. In the following analysis, we will focus on the radial range from 6 to 12 kpc for both the Galactocentric radius (R) and the guiding radius (R g ). In order to enhance the visualization of the phase space structures, at 6 and 7 kpc, the radial annulus width is set to 0.6 kpc (i.e., ±0.3 kpc), while the annulus width is 0.4 kpc (i.e., ±0.2 kpc) at the other radii. The typical velocity uncertainty is about 1 km/s for the radial, azimuthal, and vertical velocities (Gaia Collaboration et al. 2018b;Antoja et al. 2018). The guiding radius (R g ) of each star is calculated according to the rotation curve in Huang et al. (2016). To highlight the snail shell feature, we adopt the number density contrast map used in Laporte et al. (2019) and Li & Shen (2020). To evaluate the influence of the parallax bias in the Gaia catalog, we also test the same analysis with the parallax corrected Gaia sample , which contains 6,565,715 stars with relatively smaller spatial coverage. The main results in this work are robust and not affected.
3. PHASE SPACES ACROSS THE GALACTIC DISK 5 The main results are not affected if we choose other values of the solar peculiar motion, e.g., Tian et al. (2015) or Huang et al. (2015). 6 φ is the angle in the Galactocentric polar coordinate which increases clockwise. The Sun−Galactic center line has φ at 0 • .
In this section, we explore several different phase spaces of the sample, namely, the R − V φ , R − R g , V R − V φ , and the Z − V Z phase spaces. As we will show later, major structures in these phase spaces are consistent with previous studies, confirming the robustness of the sample to trace the kinematic status of the Milky Way disk.
3.1. R − V φ and R − R g Phase Spaces Fig. 4 shows the number density maps of the main sample in the R − V φ (top) and R − R g (L Z ) phase spaces (bottom). The top panel reveals several diagonal ridge lines, consistent with the previous observations Laporte et al. 2019;Fragkoudi et al. 2019). In the bottom panel, several nearly horizontal stripes can be seen across the phase space, suggesting a possible connection between the ridges in the R − V φ phase space and the constant angular momentum curves. However, as shown in Khanna et al. (2019) and Fragkoudi et al. (2019), at higher V φ the stripes may be more consistent with the constant energy lines. The observed diagonal ridges could have quite complicated origins, considering the presence of the ridge line structure in the V R or metallicity color-coded R − V φ phase space (Fragkoudi et al. 2019;Laporte et al. 2019;Liang et al. 2019;Wang et al. 2020).
From the bottom panel, stars at certain guiding radius are actually located within a large radial range. At the radius far from the solar neighborhood, the number of stars can be improved by selecting stars according to their guiding radii. This could help to enhance the phase space snail shell feature at a larger radial range. Fig. 5 shows the V R −V φ phase space of stars at different radius, with the number density and the J R color-coded maps in the top and bottom rows, respectively. From the main sample, we study the velocity phase space distributions at the 7 radial annuli centered at R = 6, 7, 8, 9, 10, 11, and 12 kpc with the same width as mentioned before. Gaia DR2 has revealed the new arch-like structures in the V R − V φ phase space for stars in the solar neighborhood, enclosing the classical moving groups (Gaia Collaboration et al. 2018a;Antoja et al. 2018;Li & Shen 2020). In fact, similar arch-like features can be seen at R = 8, 9, and 10 kpc. As the radius increases, the arches seem to systematically shift towards lower V φ , consistent with previous observational studies (Antoja et al. 2014;Ramos et al. 2018).
According to the epicycle theory, at a given radius, the azimuthal velocity difference between a star and the local circular velocity can be estimated as where κ is the epicycle (radial oscillation) frequency, x is the radial displacement from the guiding radius, and γ = 2Ω g /κ. Under the flat rotation curve assumption, γ equals to √ 2. The ratio between the maximum radial velocity and the maximum azimuthal velocity difference can be estimated, which is simply 1/γ. Adopting the best-fit potential of Model I from Irrgang et al. (2013), we use AGAMA package (Vasiliev 2019) to derive the epicycle frequencies and the action values. In Fig. 5, from left to right panels (6 to 12 kpc), the γ values are 1.32, 1.35, 1.39, 1.40, 1.41, 1.43, and 1.44, respectively. In the bottom panels for the V R −V φ phase spaces color-coded with J R , the black ellipses with this axis ratio are overplotted compared to the white curves of the constant J R contours. The white contours and the black ellipses agree well with each other.
FIG. 5.-V R −V φ phase space distributions of stars in different Galactocentric radial annulus from R = 6 to 12 kpc. The top and bottom rows show the number density and the J R color-coded maps, respectively. In the top row, arches can be seen in different radius, especially at R = 8, 9, and 10 kpc. The number of stars in each radial range is also shown in the panel. In the bottom row, the white curves represent the constant J R contours, while the black ellipses are the analytical estimation of the constant J R contours based on the epicycle theory. Similarly, we choose 7 different guiding radial annuli from 6 to 12 kpc. The V R − V φ phase space in different R g ranges are shown in Fig. 6 with the same layout as Fig. 5. Compared to Fig. 5, the distribution is relatively smooth without clear arches 7 . Similar to Fig. 5, the azimuthal velocity difference between a star and the circular velocity at its guiding radius can be estimated as V φ (R 0 ) − V circ (R g ) = γκx/2. 8 Therefore, the ratio between the maximum radial and azimuthal velocity difference is simply γ/2. The γ values from R g = 6 to 12 kpc are 1.35, 1.37, 1.39, 1.41, 1.43, 1.45, and 1.46, respectively. The axis ratio is slightly smaller than that in Fig. 5. The black ellipses with this axis ratio are shown in the bottom panels, with good agreement with the white contours.

Z − V Z Phase Space
In Fig. 7, the Z −V Z phase spaces of stars at different Galactocentric radii (R) are shown in the left figure, with the first, second and third columns representing the number density contrast (∆N), V R and V φ color-coded maps, respectively. The number density contrast is derived as ∆N = N/ N − 1, where N 7 An elongated overdensity can be recognized at Rg = 6 kpc with V φ ≈ 170 km/s, and at Rg = 7 kpc with V φ ≈ 200 km/s. These stars mainly come from the solar radius.
is the Gaussian kernel convolved number density distribution 9 (Laporte et al. 2019;Li & Shen 2020). The snail shell can be recognized in all the panels, confirming the snail shell feature as a global phenomenon across the Galactic disk. The shapes of the snail shells as revealed by the ∆N map vary at different radius. As the radius increases, the snail shell becomes more elongated along the Z axis, more squashed along the V Z axis, and less wound, due to the decreasing vertical restoring force and oscillation frequency at larger radius. These results are consistent with Laporte et al. (2019) and Wang et al. (2019).
The snail shells at 6, 7, and 12 kpc are relatively weak and hard to discern in ∆N maps, although at R = 7 kpc, the V φ color-coded phase space still reveals a snail shell pattern. At R = 8 kpc, the snail shell in ∆N map looks roughly similar to those in the V φ color-coded map, which is not true for R = 9 and 10 kpc. As already demonstrated by Li & Shen (2020), the snail shell appeared in the V φ color-coded phase space may not truly reflect the real phase mixing signal due to the variation of the snail shell shape at different V φ (also see Laporte et al. 2019).
The R g -based Z − V Z phase space maps are shown in the right figure of Fig. 7. Clear snail shells can be seen in the and guiding radial annulus (right figure). In the left (right) figure, the first, second, and third columns show the number density contrast ∆N map, V R , and V φ color-coded phase spaces, respectively. Each pixel corresponds to 0.02 kpc × 1.2 km/s. At each radius, the snail shells are more clearly visualized in the ∆N map in the guiding radial annulus in the right figure. Vertical stripes due to the observational bias become much less noticeable in the Rg-based snail shell.
∆N maps, which look similar to the snail shells in the corresponding R-based phase spaces with improved clarity. On the other hand, there is no clear snail shell-like pattern in the V R or V φ color-coded phase space shown in the second and third columns. Patterns in the V φ color-coded phase spaces are discussed more in Appendix A. In Appendix B, we show R g -based phase spaces for 30 annuli, which are evenly sampled between 6 and 12 kpc separated by 0.2 kpc each. A gradual variation of the snail shell shape towards larger radius is quite clear. In addition, in the R g -based phase space, the gap in the R-based phase space around Z = 0 (lack of stars due to dust extinction) is also mitigated. More details are shown in Appendix C.

ANGULAR MOMENTUM: A FUNDAMENTAL
PARAMETER TO TRACE THE SNAIL SHELL In the epicycle approximation, a star revolves around the guiding center which is on a circular orbit around the Galac-FIG. 8.-Phase space snail shells shapes measured at different R (top row) and Rg (middle row) from 8 to 11 kpc. The overlaid green curves with error bars represent the extracted profiles of the snail shells. Apparently, compared to the top row, more wraps of the shell can be identified in the middle row at each Rg. The bottom row compares the phase space snail shell shapes between R (blue) and Rg (red) from 8 to 11 kpc. φ and Rnorm are the angle (increasing counterclockwise in the phase space) and normalized dimensionless radius of the measured snail shell profiles in the polar coordinates. The snail shell shapes at the same R and Rg are consistent, except that the profiles in Rg extend further with larger φ which indicates more wraps of the snail shell in the phase space. tic center. Grouping stars according to the guiding radius instead of the present Galactocentric radius can avoid the mixture of stars with different angular momentum. In this section, we perform quantitative comparisons of the phase space snail shells between different R and R g ranges, and compare to other previous works to emphasize the importance of guiding radius to trace the vertical phase mixing signal.

Snail Shell Shape Measurement
As shown in the ∆N maps of the R-and R g -based phase spaces in Fig. 7, the snail shell is reflected with dark-brown color (local maximum) separated by light-brown color (local minimum). Ideally, if we work in the polar coordinate of the phase space, and to choose a wedge centered at (0, 0), then averaging the signal at each phase space radius 10 in the wedge will result in a ∆N profile, with the peaks and troughs corresponding to the snail shell and inter shell regions, respectively. After identifying the shell positions in each wedge, the 10 In the phase space, Z and V Z values at each point are divided by 1.5 kpc and 70 km/s, respectively, to get a normalized dimensionless radius Rnorm = Z 2 norm +V 2 Z,norm .
identified peak positions in all the wedges can be simply connected to recover the whole snail shell structure. The phase space is equally divided into 12 wedges with azimuthal angle of 30 • . The phase space angle φ = 0 • is defined in the direction of Z = 0 and V Z < 0 that increases counterclockwise following the snail shell outwards from the central region of the phase space. After deriving ∆N − R norm profile, the local peak and trough positions are identified; the snail shell shape is retrieved by connecting all the peak positions in the wedges. We have identified the snail shell structure in the ∆N maps at R, R g = 8, 9, 10, and 11 kpc, which are shown in Fig. 8 with the measured snail shell shape overlaid on the ∆N map. 11 The outermost wrap of the snail shell is truncated depending on its relative amplitude. At 8 and 9 kpc, the snail shell is measured up to φ = 720 • , then the measurement is terminated when the snail shell relative amplitude 12 is below 0.025. For FIG. 9.-The R-and Rg-based phase space maps at the solar radius (R, Rg = 8.34 kpc) using the Gaia RVS and the main sample (similar to Fig. 7). The R-and Rg-based phase space snail shells look similar with improved clarity and one more wrap seen in the middle and bottom rows in Rg-based phase space. 10 and 11 kpc, the snail shell is traced to φ = 540 • , then the profile is truncated when the snail relative amplitude is lower than 0.1. The error bar of the snail shell shape is determined by 1/4 separation between the two troughs on each side of the local peak (except for the outer most shell, which is determined by 1/2 separation between the peak and the inner trough since no outer trough can be determined). As shown in Fig. 8, the green curves with error bars generally follow the visually identified snail shell. Due to the large fluctuations in the signal, the outer shell is relatively difficult to map compared to the inner shell. Clearly, the green curve marking the R g -based snail shell shows more wraps and well defined snail shell shape than the R-based snail shell.
The bottom row of Fig. 8 shows the shape of the R-and R gbased snail shell in the R norm − φ plane. At each radius, the two types of points follow the same track, suggesting that the snail shell shapes are similar. The red dot (R g ) extends to larger φ and R norm values to show more wraps in the phase space. In addition, the slope in R norm − φ plane is steeper at larger radius compared to smaller ones. This is consistent with the theoretical expectation that snail shell at larger radius is more loosely wound than the snail shell at the smaller radius.

Comparison with Previous Works
We first compare our results with Antoja et al. (2018), the discovery paper of the phase space snail shell. As shown in the top row of Fig. 9, we reproduce the phase space snail shell with stars near the solar radius in the Gaia RVS sample (R = 8.34 ± 0.1 kpc), similar to Fig. 1 in Antoja et al. (2018). For comparison, in the middle row, we also generate the R g -based snail shell using Gaia RVS data. The bottom row shows the R g -based snail shell using the LAMOST + Gaia sample. As shown in the ∆N maps (second column), it seems that the snail shell features are similar except R g -based snail shells are clear with one more wrap.
The snail shell shapes of the two phase spaces are measured with the method mentioned in the previous section. The ∆N maps with the snail shell curves overlaid are shown in the left three panels in Fig. 10. The comparison between the three curves are shown in the right panel. Clearly, the snail shell shapes are consistent within φ < 700 • (in the inner region of the phase space). The red and green dots (R g ) extend to larger R norm and φ values.
The snail shell clarity is also estimated with the ∆N profiles extracted along Z and V Z stripes; in the ∆N map, we FIG. 12.-The Z −V Z phase space distributions of stars in different angular momentum ranges. The median L Z increases from 1500 to 2400 kpc km/s from left to right columns, with the number density, density contrast, and L Z color-coded phase spaces from top to bottom rows. In each panel, the bin width ∆L Z is 200 kpc km/s (±0.4 kpc in Rg) to be consistent with Khanna et al. (2019). Panels in the bottom row well reproduce Fig. 16 in Khanna et al. (2019). The difference between the L Z color-coded phase space and the ∆N map is quite clear. choose the stripe with |V Z | < 15 km/s to get the ∆N − Z profile, and the stripe with |Z| < 0.15 kpc to get the ∆N −V Z profile. The results are shown in Fig. 11 with the left and right panels corresponding to the ∆N profiles along Z and V Z , respectively 13 . The yellow vertical solid and dashed lines mark the snail shell turning-around points (V Z = 0) and mid-plane points (Z = 0), respectively. The local peak position in the red and green curves (R g ) agree quite well with the yellow lines, with relatively larger differences between the peak and 13 The profiles are truncated with ∆N uncertainty larger than 0.14. trough positions (∼ 0.08). The blue curve (R) also seems to agree with the position of the snail shells, but with relatively smaller difference between the peaks and troughs (∼ 0.05). The comparison confirms that in the solar neighborhood, the snail shell in the R g -based phase space is better revealed than the traditional R-based phase space.
Based on a small sample covering ∼ 1 kpc distance from the Sun, Khanna et al. (2019) found that the phase spiral pattern for a given L Z is almost invariant with radius, suggesting the angular momentum as a more robust quantity to characterize the snail shell compared to V φ . As shown in Fig. 16  FIG. 15.-Ω Z − √ J Z distributions of stars at R = 8 kpc (top row) and Rg = 8 kpc (middle row). The top left and right panels show the number density and V φ color-coded Ω Z − √ J Z distributions at R = 8 kpc, respectively. The middle left and right panels show the number density and ∆Rg = Rg − 8 kpc color-coded Ω Z − √ J Z distributions, respectively. The distribution is much narrower for stars in different guiding radius bin (middle row), resulting in the clear phase space snail shell. At certain radius R, Ω Z generally decreases with increasing V φ , producing the snail shell in the V φ color-coded phase space. At Rg = 8 kpc, Ω Z decreases slightly with increasing ∆Rg. The bottom left panel shows the Ω Z − √ J Z distribution of stars at Rg = 8 ± 0.01 kpc with J R < 0.01 (blue) and J R > 0.06 (red). The bottom right panel compares the Ω Z − √ J Z distributions of stars at different Rg (from 6 to 11 kpc with 0.02 kpc width annulus) with J R < 0.01. Stars that closely follow circular motions at each guiding radius (with small J R ) show very tight correlations between Ω Z and √ J Z .
of Khanna et al. (2019), the snail shell pattern revealed in the L Z color-coded phase space (∆L Z = 200 kpc km/s) with the orientation of the snail shell changing with L Z . Our result improves upon Khanna et al. (2019) in the aspect that we have explored a large radial range and performed comparisons of the ∆N phase spaces (representing the intrinsic shape of the snail shell rather than V φ or L Z color-coded phase space) of different R and R g ranges. According to our results, R g (or L Z ) is not only robust, but fundamental to trace the phase space snail shell across the Galactic disk.
To compare our results with Khanna et al. (2019), we explore the phase space distributions at each L Z . We first derive the 6D information for stars in the Gaia RVS sample using the same solar position and kinematic configuration as Khanna et al. (2019). Then we are able to reproduce Fig. 16 in Khanna et al. (2019) by selecting stars in different angular momentum bin with ∆L Z = 200 kpc km/s (∼ ±0.4 kpc in R g ). Fig. 12 shows the number density, ∆N, and L Z color-coded phase spaces from top to bottom rows, respectively. Apparently, the snail shell shapes revealed in the top and middle rows (N and ∆N maps) are quite different from the corresponding L Z color-coded maps in the bottom row. This confirms our previous conclusion that it is the number density map representing the intrinsic phase space structures. The L Z color-coded phase space may not accurately reflect the phase space snail shell structure because it is the number weighted average of the angular momentum of stars. In an imaginary extreme case, where all the stars in each panel of the top row in Fig. 12 have the same angular momentum, the L Z colorcoded Z − V Z phase space will not reveal any feature at all. A good example can be seen in Fig. 25 in Appendix B, where ∆L Z = 100 kpc km/s (much smaller than the value adopted by Khanna et al. 2019) and no snail shell-like pattern even exists in the L Z color-coded phase space at different radius. We can conclude that R g -based snail shells can reveal more clear, and intrinsic information on the phase mixing process than the L Z color-coded ones.

Colder and Hotter Orbits Dichotomy
Li & Shen (2020) found that, for stars near the solar neighborhood, clear phase space snail shell can only be seen in the dynamically colder orbits (|V φ − V LSR | < 30 km/s or J R < 0.04), while the hotter orbits have probably phase-wrapped away already. The perturbation should happen at least ∼ 500 Myr ago to facilitate the phase mixing of the hotter orbits.
Here with a larger sample we also investigate the cold/hot dichotomy across the Galactic disk. Similar to Li & Shen (2020), at different radius, we use the same criteria to separate the colder and hotter orbits. The phase space distributions for the colder and hotter orbits at different radius are shown in the  FIG. 17.-The evolution of the test particle simulation with the impulse approximation implemented. The intruder is assumed to hit the disk at 18 kpc along the X−axis from the center, with a vertical downward velocity of 300 km/s. From left to right, the four columns show the evolution of the model at T = 50, 250, 500, and 800 Myr after the impact. We adopt the snapshot at T = 500 Myr in the phase space analysis.
left and right figures of Fig. 13, respectively. Consistent with Li & Shen (2020), prominent snail shell can only be seen on the colder orbits at different radius.
We also test to divide stars at each guiding radius into colder and hotter orbits. The Z − V Z phase space distributions are shown in Fig. 14 with the colder and hotter orbits shown in the left and right figures, respectively. In fact, we have raised the orbital hotness criteria to J R > 0.06, since snail shell features are still visible in the warmer orbit at intermediate J R values (0.04 < J R < 0.06). This test again confirms the importance of the guiding radius in revealing the phase space snail shell features, with the radial orbit hotness (ellipticity) playing a less important role.

Vertical Action and Oscillation Frequency Analysis
Aiming to better understand the vertical phase mixing process, we explore the action space of the sample. Fig. 15 shows the Ω Z − √ J Z correlations for stars at R, R g = 8 kpc. At each radius, the anti-correlation between Ω Z and √ J Z is the main reason behind the phase space snail shell, representing the anharmonic vertical oscillation. Comparing the top and middle panels, the Ω Z − √ J Z distribution at each R g is much tighter than that at different R, resulting in a clear snail shell in the number density map of R g selected sample. This also implies that regardless of the current position in the disk, stars with the same angular momentum tend to follow similar Ω Z − √ J Z correlation (and similar snail shell shape in the Z − V Z phase space).
The top right panel in Fig. 15 shows the same Ω Z − √ J Z correlation but color-coded with V φ at R = 8 kpc. A clear trend can be seen, that Ω Z decreases for larger V φ at a given J Z , consistent with Binney & Schönrich (2018). On the other hand, the middle right panel shows the Ω Z − √ J Z correlation of stars at R g = 8 kpc color-coded with the guiding radius deviation (∆R g = R g − 8 kpc). The relatively smaller guiding radius corresponds to higher vertical oscillation frequency Ω Z .
To better understand the tight correlation between Ω Z and √ J Z in the R g -based sample, we first choose stars in a very narrow guiding radius range, i.e., R g = 8 ± 0.01 kpc, which are further separated into a dynamically colder subsample with J R < 0.01 and a dynamically hotter subsample with J R > 0.06.
The Ω Z and √ J Z distributions of the two subsamples are shown in the bottom left panel of Fig. 15. The correlation is very tight for the dynamically colder orbits at very small J R , with the dispersion increasing considerably at larger J R for the dynamically hotter orbits. This is expected since the dynamically colder stars are close to circular; these stars with different vertical action move in the vertical gravitational potential well at the same radius, thus resulting in tight correlation between Ω Z and J Z . On the other hand, the dynamically hotter stars have large radial excursion, with the corresponding vertical profile of the gravitational potential well varying to result in the large dispersion. Notice that the dispersion is mainly towards lower vertical frequency. This is probably due to the fact that the dynamically hotter stars spend more time around apocenter in the outer disk, with the vertical oscillation frequency reduced in the shallower potential well. The bottom right panel in Fig. 15 shows the tight Ω Z − √ J Z correlation for dynamically colder stars at different guiding radius from 6 to 11 kpc, each with 0.02 kpc width annulus and J R < 0.01. At a given vertical action value, the vertical oscillation frequency decreases progressively at larger guiding radius. The slope at larger R g is also shallower, consistent with the loosely wound snail shell shown in Fig. 7.

TEST PARTICLE SIMULATIONS
To better understand the spatial variation and the importance of the angular momentum of the phase space snail shell across the Galactic disk, we perform test particle simulations with 40,000 test particles rotating in a realistic Milky Way potential (Irrgang et al. 2013). Initially, the test particles are evenly distributed in the disk mid-plane from R = 6 to 12 kpc with no vertical motion (Z = 0,V Z = 0). The distribution of the azimuthal velocity follows the Gaussian distribution that peaks at the local circular velocity with dispersion of 10 km/s. The median radial velocity is set to 0 with dispersion of 10 km/s. These test particles are evolved for 10 Gyr to be fully phase mixed in the horizontal direction. Afterwards, we impose the vertical perturbation on all the test particles in two different approaches. The first approach is similar to Li & Shen (2020), where the particles receive vertical downward velocities with the vertical positions barely changed (median Z = 0 and V Z = −20 km/s with a dispersion of 0.02 kpc and 10 km/s, respectively). For the second ap- FIG. 18.-The Z −V Z phase space distributions of stars in different radial ranges (from R = 6 to 12 kpc) of the test particle simulation in the second approach. The left three columns show the number density maps in the Z − V Z phase space at each radius (different rows) and each azimuthal wedge (different columns). The middle three columns show the V R color-coded the Z −V Z phase spaces, and the right three columns show the V φ color-coded phase spaces corresponding to these number density maps. proach, we consider a more realistic impulse approximation on the test particle velocities following Binney & Schönrich (2018), to mimic the effect of a fly-by external perturber.
In the first approach, regardless of the different azimuthal angles, the test particles at the same radius receive the same vertical perturbation, and thus the same vertical velocity distributions. In the following analysis, all the particles in each radial (or guiding radial) annulus are used. After 500 Myr evolution, the Z − V Z phase space distributions are shown in Fig. 16 for different radial annuli (left figure) and different guiding radial annuli (right figure) similar to Fig. 7. The Rbased snail shell shapes are similar to R g -based snail shell, with the R g -based snail shells narrower and well-defined than the R-based ones. These results agree with the observational results in Fig. 7. Stars with the same angular momentum tend to follow the same vertical oscillation pattern (although they were located at different radius). On the other hand, in a given radial annulus, stars usually have a wide distribution of the angular momentum. The mixture of stars with different angular momentum (and the related phase space pattern) results in the blurring of the snail shell.
In the left figure of Fig. 16, the third column shows the V φ color-coded phase space. Different color shows different snail shell patterns. For particles at a given radius, the azimuthal velocity is proportional to the angular momentum.
Since stars with different angular momentum have different snail shell pattern, when color-coded with V φ , different snail shells naturally emerge in the phase space. Comparing the results with different R g ranges in the right figure, there is no systematic variation of the snail shell shape with azimuthal velocity, since the median V φ at different part of the phase space along the snail shell is close to 0. In Fig. 16 for both the R-and R g -based phase spaces, no variation of the snail shell pattern with V R can be seen. The median radial velocity in the phase space is close to 0. This is mainly due to the initially imposed symmetric distribution of V R in our test particle simulation.
For the second approach, the perturbation to the vertical velocities of the particles are estimated following the impulse approximation in Binney & Schönrich (2018). The impact parameter of the external perturber is 18 kpc along the X−axis from the Galactic center at a speed of ∼ 300 km/s. The mass of the perturber is 2 × 10 10 M ⊙ and T = 66 Myr as the passage timescale. To work in the non-inertial frame with the Galactic center stationary, the gravitational pull from the intruder to a point mass at the Galactic center is subtracted from the acceleration of the intruder on the other test particles. The in-plane velocity change δv is computed by multiplying the characteristic timescale of the passage with FIG. 19.-The Z −V Z phase space distributions of stars in different guiding radial ranges (from Rg = 6 to 12 kpc) of the test particle simulation in the second approach. The left three columns show the number density maps in the Z − V Z phase space at each guiding radius (different rows) and each azimuthal wedge (different columns). The middle three columns show the V R color-coded the Z − V Z phase spaces, and the right three columns show the V φ color-coded phase spaces corresponding to these number density maps. the acceleration. A downward component δv ⊥ is added as δv ⊥ = αδv (1 − β sin (φ − φ intruder )) with α = 0.4 and β = 0.5 (see Eqt. 3 in Binney & Schönrich 2018). The vertical perturbation is composed of two components: a constant term associated with the standard dynamical friction and a term proportional to sin (φ − φ intruder ). The second term arises from the consideration that stars initially moving towards the intruder (φ − φ intruder < 0) receive a net downward kick, while those stars initially moving away from the intruder (φ − φ intruder > 0) receive a net upward kick. We define φ intruder = 0 • along the X− axis, and the azimuthal angle increases counterclockwise.
After the impact, the time evolution of the test particles is shown in Fig. 17. As shown in the first column, 50 Myr after the perturbation, the disk becomes lopsided that is stretched towards the intruder with a prominent single arm. The amplitude of the vertical perturbation of each particle depends on the azimuthal angle, with larger amplitude for particles closer to the intruder. This is the main difference compared to the previous test particle simulation where the imposed vertical perturbation is azimuthally symmetric. In the second column, 250 Myr after the impact, the initial perturbation gradually evolves to form the m = 2 mode due to the differential rotation, with the vertical oscillation in different parts of the disk. After 500 Myr evolution (third column), the m = 2 mode is still present, with continuous phase mixing in the vertical di-rection of the disk. In fact, this is the snapshot we adopt in the following analysis, since the phase space snail shell is most prominent and roughly similar to the observations. At later times, as shown in the fourth column (800 Myr), the face-on image still exhibit prominent m = 2 mode feature, with the roughly symmetric distribution of the particles with respect to the disk mid-plane in the edge-on view of the model.
Recently, Bland-Hawthorn & Tepper-Garcia (2020) used a high resolution N-body model to simulate the impact of an external intruder on the main disk. They found good agreement with the model in Binney & Schönrich (2018) indicating that the impulse approximation is valid to describe such process. After the impact, Bland-Hawthorn & Tepper-Garcia (2020) noticed a strong, m = 1 bending mode across the disk is set up, which gradually wraps up due to the differential rotation of the disk to result in a superposition of two distinct bisymmetric (m = 2) modes, namely a spiral pattern and a bending wave. Their self-consistent N-body simulation results are quite similar to our simple test particle simulation under the impulse approximation.
Since the vertical perturbation is not azimuthally symmetric, different phase space snail shells are expected to emerge at different azimuthal angles. The disk is separated into three azimuthal wedges with 120 • each, in order to trace the azimuthal variation of the phase space snail shell with sufficient number of particles. The Z − V Z phase space distributions at different radius and azimuthal wedges are shown in Fig. 18. At each radius, as expected, the shape of the phase space snail shell is different in different azimuthal wedge. For example, at R = 8 and 9 kpc, the snail shell is barely seen in the azimuthal wedge with φ = (4π/3, 2π) (first column), while in the other azimuthal wedges at (2π/3, 4π/3) (second column) and (0, 2π/3) (third column), the snail shells are clear. In the V R (middle three columns) and V φ (right three columns) colorcoded phase spaces, evidence of snail shell shape changing at different velocities can be seen, roughly consistent with observations.
The phase space distributions at different guiding radius are shown in Fig. 19. With azimuthal variations still present, the R g -based number density phase space (left three columns) show prominent snail shells compared to the R-based phase spaces in Fig. 18. In the V R and V φ color-coded phase spaces, there is no systematic variation of the snail shell shape with the velocities.
Bland-Hawthorn & Tepper-Garcia (2020) failed to detect the phase space snail shell feature in the angular momentum (guiding radius) bins, indicating that the N-body model is still incomplete to match all the observational results. More efforts are needed in both observations and simulations to better understand the vertical phase mixing process in the future.
6. SUMMARY We utilize ∼ 11 million stars from both the Gaia RVS sample and the LAMOST DR6 sample to study the vertical phase mixing across the Galactic disk. We confirm the existence of the snail shells in the Z − V Z phase space in the Galactic disk from 6 to 12 kpc, and quantitatively measure the shape of the snail shells. We find that the guiding radius (angular momentum) is fundamental to reveal the vertical phase mixing signals. Compared to the R-based phase spaces, the clarity of the snail shell in the R g -based phase space is increased; more wraps of the snail shell can be seen; phase space is less affected by the lack of stars close to the disk mid-plane due to extinction; snail shell signal is also amplified in a greater radial range.
Additional evidence comes from the cold/hot dichotomy discovered in Li & Shen (2020) in the solar neighborhood. For the R-based phase space, the cold/hot dichotomy still holds, with prominent snail shell only in the colder orbits. However, in each R g range, snail shell can still be seen in warmer orbits. The snail shell eventually disappears for very hot orbits with J R > 0.06 (compared to the original criteria J R > 0.04 for the hotter orbits in Li & Shen 2020). Therefore, to study the vertical phase mixing process, it is recommended to use R g (with additional constraints on J R ) to better reveal the intrinsic shape of the snail shell in the Galactic disk. In another ongoing work, we are trying to model the snail shell shapes at different radius to constrain the vertical shape of the Galactic potential.
The difference between the phase space distributions in the Galactocentric radial range and the guiding radius range can be understood in the Ω Z − √ J Z plane. At each radius, the distribution of stars in the Ω Z − √ J Z plane shows anticorrelation with broad distribution, while stars in the corresponding R g range show a much tighter correlation, resulting in more prominent snail shells in the phase space.
Test particle simulations are also performed to understand this phenomenon. In the first approach, a simple azimuthally symmetric vertical perturbation is imposed on all the test particles. Snail shell features can be seen across the disk, which are more prominent when the particles are grouped into different guiding radius ranges. In the second approach, the perturbation with the impulse approximation is applied on all the test particles. Again, the phase space snail shell is more prominent in different guiding radial ranges, consistent with the observations. I thank Juntai Shen for insightful suggestions that help improve the clarity of the paper, and Chao Liu for helpful discussions. The research presented here is supported by the "111" Project of the Ministry of Education under grant No. B20019, and by the MOE Key Lab for Particle Physics, Astrophysics and Cosmology. This work made use of the Gravity Supercomputer at the Department of Astronomy, Shanghai Jiao Tong University, and the facilities of the Center for High Performance Computing at Shanghai Astronomical Observatory.
Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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.

APPENDIX
A. V φ COLOR-CODED PHASE SPACE PROPERTIES As shown in the right figure in Fig. 7, in the V φ color-coded phase space at R g = 8, 9, 10, 11, and (possibly) 12 kpc, the average azimuthal velocity seems smaller at larger |Z| and higher at larger |V Z |. This feature can be understood from the statistical distribution of stellar orbital properties. Focusing on R g = 8, 9, 10, and 11 kpc, we select three subsamples at each guiding radius bin, with subsample 1 at |Z| ≤ 0.5 kpc and |V Z | ≤ 15 km/s, subsample 2 at |Z| ≤ 0.5 kpc and |V Z | ≥ 20 km/s, and subsample 3 at |Z| > 0.7 kpc. Fig. 20 shows the distributions of J Z , Z, R and V φ at each guiding radius, with the black solid line, the red dotted line and the blue dashed line representing the subsamples 1, 2, and 3, respectively. As expected from the location in the Z −V Z diagram, and also shown in the first column, the vertical action is smallest in subsample 1, and largest in subsample 3. The vertical height is largest in subsample 3 (second column), corresponding to larger radial range shown in the third column; stars with larger vertical hights tend to be located outwards at larger radius where the disk potential well is shallower. At each guiding radius, these three subsamples have similar angular momentum. Therefore, the larger radii naturally result in smaller azimuthal velocities in subsample 3 (fourth column).
For subsample 2, although the vertical action is higher than that of subsample 1 as shown in the first column, they show similar vertical height distributions (second column), with subsample 1 slightly narrower. The much larger vertical velocities of subsample 2 hint for smaller radius for this subsample; stars are more likely to acquire larger vertical velocities when traveling at smaller radii with deeper potential well. This picture is consistent with the third column. Therefore, the azimuthal velocity of subsample 2 is relatively higher than that of subsample 1.
Besides the V φ pattern of the phase space at each R g , careful inspection of the left columns in Fig. 7 for the snail shell at different radial annuli indicates differences between the snail shell shapes of the ∆N map (first column) and V φ color-coded phase spaces (third column), especially for the annulus at 10 kpc. This difference is likely caused by the snail shell shape variation at different V φ . According to Li & Shen (2020), from the mathematical point of view, the V φ color-coded phase space can be considered as the number weighted average of the azimuthal velocity of stars in each radius. With V φ color-coding, problems with the incomplete sampling close to the mid-plane can be alleviated. However, the V φ color-coded phase space may not truly reflect the real phase space snail shell shape. Li & Shen (2020) has shown in their Fig. 14 that by combining stars following different snail shell shapes at slightly different V φ , the phase space color-coded with V φ would look quite different from the ∆N map. As shown in the second and third columns of the right figure for different guiding radius ranges, there is even no well recognized snail shells at all for the V R and V φ color-coded phase spaces.

B. Z − V Z PHASE SPACE AT DIFFERENT R G
In order to better visualize the gradual transition of the phase space snail shell across the Galactic plane, we separate the sample into 30 guiding radius bins from 6 to 12 kpc with the bin width of 0.2 kpc. The number density map and the contrast maps are shown in Figs. 21 and 22, respectively. As the guiding radius increases, the snail shell becomes less wound and more extended along the Z direction. The gradual transition is very clear. The V R , V φ , and L Z color-coded phase spaces are shown in Figs. 23, 24, and 25, respectively. No clear pattern can be seen, consistent with our previous results.

C. THE GAP IN THE PHASE SPACE AROUND Z = 0
In Fig. 7, a clear gap close to Z = 0 in the R-based phase space number density map can be seen, especially at radius smaller than 8 kpc. This feature is mainly caused by the stronger dust attenuation in the disk mid-plane. However, this gap becomes much weaker or even disappears in the corresponding phase space at R g . Here we perform detailed comparison between the LAMOST and Gaia RVS sample at each R and R g to better understand this difference. Fig. 26 shows the comparison at 8 kpc. The gap is very clear in the LAMOST sample at R = 8 kpc (top row), and much weaker at R g = 8 kpc (middle row). For the Gaia RVS sample, the gap is weaker, but the improvement in R g is apparent. The bottom panel shows the same R g phase space as the middle row but color-coded with ∆R = R − 8 kpc. It seems that at R g = 8 kpc, most stars come from R = 8 ∼ 9 kpc. In this case, stars in the disk are reshuffled with a large fraction of stars near the solar radius (where the observed stellar number density is highest close