Measuring the Milky Way Vertical Potential with the Phase Snail in a Model-independent Way

The vertical phase-space spiral (snail) is a direct sign of disequilibrium of the Milky Way’s disk. Nevertheless, the wrapping of the phase snail contains information on the vertical potential. We propose a novel method to measure the vertical potential utilizing the intersections between the snail and z/V z axes, for which we know the maximum vertical heights (Z max) or velocities ( Vz,max ). Using a refined linear interpolation method, we directly obtain (Zmax,12Vz,max2) for these snail intersections to constrain the vertical potential profile empirically. Our method is model-independent, since no assumptions about the snail shape or the vertical potential have been made. Although the snail binned by the guiding center radius (R g ) is more prominent, it traces a vertical potential shallower than that of the snail binned by the same galactocentric radius (R). We apply an empirical method to correct this difference. We measure the snail intersections in several R g bins within 7.5 < R g < 11.0 kpc for Gaia DR3 and apply the interpolation method to deduce the potential values at several vertical heights. The potential at the snail intersections, as well as the following mass modeling, is consistent with the popular Milky Way potentials in the literature. For the R g -binned phase snail in the solar neighborhood, the mass modeling indicates a local dark matter density of ρ dm = 0.0150 ± 0.0031 M ⊙ pc−3, consistent with previous works. Our method could be applied to larger radial ranges in future works to provide independent and stronger constraints on the Milky Way’s potential.


INTRODUCTION
A main goal of dynamical studies of our Milky Way is to investigate the mass distribution, consisting of a thin disk, thick disk, bulge/bar, stellar halo, and a dark matter halo (Bland-Hawthorn & Gerhard 2016).The different characteristics and interlinks of the various Galactic components could reveal the formation and evolution history of our Milky Way (e.g.Read et al. 2008;Helmi et al. 2018;Mackereth et al. 2019;Helmi 2020).For example, the shape of the Galactic dark matter halo helps to understand the influence of the baryonic effects to the structure formation, and thus the current galaxy formation models (e.g.Kazantzidis et al. 2010;Koposov et al. 2010;Lux et al. 2012;Eadie & Jurić 2019;Fardal et al. 2019;Iorio & Belokurov 2019;Wegg et al. 2019).In the Solar neighborhood, many dynamical studies aim to measure the local dark matter density, which is an Corresponding author: Zhao-Yu Li, Juntai Shen lizy.astro@sjtu.edu.cn;jtshen@sjtu.edu.cnimportant physical parameter to provide guidance for the direct dark matter detection experiments (e.g.Fairbairn et al. 2013;Green 2017;Maity et al. 2021;Bechtol et al. 2022).
Various methods have been applied to measure the Galactic potential.In the local region, the vertical potential and mass distribution have been estimated by analyzing the kinematics of different stellar tracers, employing the vertical Jeans equation or the distribution function method (e.g.Zhang et al. 2013;Xia et al. 2016;Hagen & Helmi 2018;Schutz et al. 2018;Sivertsson et al. 2018;Buch et al. 2019;Guo et al. 2020;Salomon et al. 2020).The local baryonic census studies help to clarify the baryonic contribution to the Galactic potential and provide appropriate priors to the modeling (e.g.Holmberg & Flynn 2000;Flynn et al. 2006;McKee et al. 2015).Some of more global studies of the Galactic mass distribution usually fit a specific galactic mass model to the Galactic rotation curve, which is obtained by the compilation of kinematic measurements of gas, stars and masers (e.g.Pato et al. 2015;Huang et al. 2016;Benito et al. 2019Benito et al. , 2021;;de Salas et al. 2019;Karukes et al. 2019;Lin & Li 2019;Sofue 2020), or by the radial Jeans equation (e.g.Eilers et al. 2019).According to the recent measurements, the local dark matter density is estimated within 0.4 − 0.6 Gev cm −3 , corresponding to 0.010 − 0.015 M ⊙ pc −3 , with uncertainties mainly contributed by the observations and the degeneracy between the baryon and dark matter density profiles (Read 2014;de Salas & Widmark 2021).
In other global mass profile studies of Milky Way, tracers such as halo stars, stellar streams, globular clusters and satellites, have also been extensively utilized, especially to study the dark matter halo shape (e.g.Fritz et al. 2018;Malhan & Ibata 2019;Hattori et al. 2021;Ibata et al. 2021;Bird et al. 2022).More recently, the action-angle based distribution function method (Binney 2010;Binney & McMillan 2011;McMillan & Binney 2013) has been employed to provide self-consistent constraints on the Milky Way potential (e.g.Bovy & Rix 2013;Piffl et al. 2014;Li & Binney 2022;Binney & Vasiliev 2023).In addition, several Milky Way potential models have been constructed through synthesizing all of the knowledge about the components of the Milky Way into a coherent picture of the gravitational potential (e.g.Irrgang et al. 2013;Bovy 2015;McMillan 2017;Cautun et al. 2020).These potentials have been widely used in the orbit integration and orbital properties estimation, in galpy (Bovy 2015) and AGAMA (Vasiliev 2019).(For reviews of local dark matter density measurements, see Read 2014 andde Salas &Widmark 2021, and for global measurements of Milky Way mass profile, see Wang et al. 2020.)Most of the dynamical studies are based on the equilibrium assumption, assuming a time invariant distribution function.However, many signs of dis-equilibrium have been found recently in Milky Way.The number density difference between the north and south of the Galactic plane displays a wavelike pattern, showing about 10% difference at some vertical heights (e.g.Widrow et al. 2012;Yanny & Gardner 2013;Wang et al. 2018a;Bennett & Bovy 2019).The disc shows a non-zero vertical bulk motion, classified into bending and breathing modes in different regions (e.g.Widrow et al. 2012;Carlin et al. 2013;Carrillo et al. 2018;Gaia Collaboration et al. 2018b;Wang et al. 2018bWang et al. , 2020;;Bennett & Bovy 2019).Recently, Antoja et al. (2018) found clear spiral features in the vertical phase space (z − V z ) color-coded by the stellar number density, azimuthal velocity and radial velocity, using Gaia DR2 (Gaia Collaboration et al. 2018a).The phase-space spiral (phase snail) is the result of incomplete phase mixing in the Galactic vertical direction.The phase snail is found to vary with stellar age, chemistry, orbital properties and the location in the Galactic disk (radially and azimuthally) (e.g.Tian et al. 2018;Bland-Hawthorn et al. 2019;Laporte et al. 2019;Wang et al. 2019;Li & Shen 2020;Xu et al. 2020;Li 2021).The phase snail has been explained by the passage of a satellite galaxy, such as the last pericentric passage of the Sagittarius dwarf (e.g.Antoja et al. 2018;Binney & Schönrich 2018;Bland-Hawthorn et al. 2019;Laporte et al. 2019;Bland-Hawthorn & Tepper-García 2021), by the buckling instability of the bar (Khoperskov et al. 2019), by the dark matter wake induced by a satellite passage (Grand et al. 2023), by cumulative effect of smaller intersections with subhalos (Tremaine et al. 2023), or by multiple origins as indicated by the two armed phase spirals in the inner Milky Way (Hunt et al. 2022;Li et al. 2023).
The observations of the north-south asymmetries and the phase snail contradict the equilibrium assumption in the dynamical modeling.Several works have tried to reconstruct the vertical potential with the consideration of the phase snail.Li & Widrow (2021, 2023) applied steady-state dynamical models to giant stars and found the phase snail in the residual of the best-fit model of the distribution function in the vertical phase space.In a series of works by Widmark et al. (2021aWidmark et al. ( ,b, 2022a,b),b), the authors predicted the snail feature by the one-dimensional (1D) phase-space orbital integration, and treated it as a component having a constant number density fraction relative to the bulk distribution function.By fitting the 2D phase space distribution in the observation, they estimated the vertical potential embedded in the phase-space orbital integration.In our previous work (Guo et al. 2022), we applied a geometric snail model superimposed on a smooth background in the phase space to fit the northsouth asymmetry of the vertical velocity dispersion profile.Our model successfully explained the north-south asymmetries in the vertical profiles of both the stellar number density and the vertical kinematics.We also obtained an estimation of the background potential which constructs the phase space distribution of the background component through the vertical Jeans equation.Conversely, under a certain Milky Way potential, the phase snail can be utilized to deduce the perturbation time through measuring the slope in the vertical frequency and action angle space (e.g.Darragh-Ford et al. 2023;Frankel et al. 2023;Widrow 2023).
The phase snail shape curve contains the information of the vertical potential.Once perturbed vertically, stars would oscillate in the vertical direction under the vertical potential, if ignoring the coupling between the vertical motion and the in-plane motion (which is a valid approximation in low vertical regions).The snail then wraps up due to the anharmonic oscillation induced by the vertical potential.In this work, we propose a model independent method to measure the Galactic vertical potential, which utilizes the intersections between the phase snail and the z/V z axes.The intersections have either known maximum vertical heights (Z max ) along the zaxis, or maximum vertical velocities (V z,max ) along the V zaxis.A proper interpolation method can be used to provide accurate (Z max , V z,max ) values for these intersections, which can be used as a direct measurement of the vertical potential values at different vertical heights.
The paper is organized as follows: The method is described and verified with a test particle simulation in Section 2. In this section, we also study the shape difference between the phase snails binned by the guiding center radius (R g ) and by the Galactocentric radius (R), and propose an empirical correction to account for this difference when measuring the vertical potential using the R g -binned phase snails, which is more prominent to be measured accurately.Section 3 describes the observational data used.The measurement of the snail intersections at different R g , the comparisons with popular Milky Way potentials and the further mass modeling utilizing the model independent potential (Z max , 1 2 V 2 z,max ) of the snail intersections are shown in Section 4. Discussions and conclusions are made in Sections 5 and 6, respectively.

Method description
The phase spiral is the result of incomplete phase mixing in the vertical phase space (z − V z space).The phase-space distribution should be steady and smooth for an equilibrium state.However, when the disc is perturbed (externally or internally), stars are disturbed to similar phase angles.As stars have different orbital frequencies, they will gradually wind up to form a spiral structure due to the differential rotation in the phase space.The wrapping of the phase snail and the dynamics behind have been described in great detail in sections 2.1 -2.3 of Widmark et al. (2021a).We briefly summarize the picture here and introduce our different method utilized for measuring the Milky Way vertical potential in this work.
The basic assumption of measuring the Milky Way vertical potential with the phase snail is that the vertical and in-plane motions are separable.More specifically, the vertical energy E z 1 is conserved for a vertically oscillating star, which moves in the vertical phase space following a certain closed trajectory with an oscillation frequency (Ω z ) and an orbital period (P(E z )).The vertical energy of a star is equal to where the vertical potential energy Ψ(z) is a reduced energy in a Milky Way potential Φ(R, z) as Given a certain value of E z , the phase-space trajectory of a star can be fully determined under certain potential model Ψ(z).The trajectory intersects with z axis at ±Z max (i.e. the turning-around points with V z = 0), and with V z axis at ±V z,max 1 Energy Ez here is in fact the energy per unit mass.
(i.e. the mid-plane points with z = 0), where the star has the maximum vertical potential and vertical kinetic energy, respectively.Thus, the vertical energy can be related to Z max and V z,max as follows: The orbital period P(E z ) of the star can be calculated by the 1D orbital integration as where the factor 4 on the right hand side is due to the symmetry about the phase space origin.The vertical frequency is therefore defined as Ω z = 2π/P (E z ).The phase-space trajectory of a star is also determined by a similar integration as Eq. 4 but with different upper limits of |z|.
As shown in the sketch plot in Fig. 1, the phase snail wraps up in the phase space and intersects with both z and V z axes.The intersections on the z axis, i.e. the odd points, reach their maximum vertical excursion with measurable Z max .Similarly, the even intersections on the V z axis have measurable V z,max , which corresponds to the vertical energy E z .Besides, the intersections have increasing vertical energy E z and decreasing vertical frequency Ω z from inner to outer regions of the phase space, with Z max and V z,max increasing monotonically.Thus, we can derive the Z max for an intersection on the V z axis by an interpolation method utilizing the two adjacent intersections on the z axis with vertical action angle difference ±π/2.For example, for point 2, its maximum vertical height can be determined by Z max,2 = f (Z max,1 , Z max,3 ), where f is an arbitrary interpolation function.Similarly, we can derive the maximum vertical velocities for intersections on the z axis with another interpolation function g, e.g. for point 3, V z,max,3 = g(V z,max,2 , V z,max,4 ).
Another important feature of the phase space motion is the linear relation between the vertical frequency Ω z and the vertical action angle θ z : where θ 0 and T evol are the initial phase angle and the evolution time since the perturbation, respectively.The adjacent intersections on the z/V z axes have a constant action angle difference of ∆θ z = π/2, which results from a constant vertical frequency difference ∆ Ω z between the adjacent intersections.Thus, a parameter having a better linear relation with Ω z would be preferred for the interpolation.For the intersections on the V z axis, both V z,max and E z can be used for the interpolation.As shown in the later Section 2.2.3, V z,max is better than E z , and a linear interpolation function would be accurate enough and model independent.After the (Z max , V z,max ) pairs have been derived for the intersection points, the relation between Z max and Ψ(Z max ) = 1 2 V 2 z,max can be built at different radii, leading to a direct measurement of the Milky Way's vertical potential.Though we have utilized the same description of the phase space wrapping as Widmark et al. (2021a) (the potential measurement of which is embedded in the orbital integration as Eq.4), we do not make any assumption of the prior potential model as they do.Our measurement would be valuable to independently assess the popular Milky Way potential models, or to constrain the Milky Way global mass model.

Method verification with test particle simulation
To verify this model independent method, we run a test particle simulation similar to Li (2021).The initial condition sampling and 3-dimensional (3D) orbital integration are performed with AGAMA (Vasiliev 2019).The input potential is the Model I from Irrgang et al. (2013).We sample a thin disc tracer population utilizing the QuasiIsothermal distribution function (DF) (Binney 2010;Binney & McMillan 2011).The basic properties such as the radial and vertical scale lengths, velocity dispersion are set to represent a dynamically cold thin disk (Binney & Piffl 2015;Bland-Hawthorn & Gerhard 2016;Li 2021).Afterwards, the disc is perturbed following the impulse approximation by imposing a vertical downward velocity with the vertical position unchanged.The vertical velocity kick received by all the particles follows a Gaussian distribution of −30 ± 5 km s −1 , i.e. the first approach used in Li (2021).The snapshot evolved for 500 Myr after the perturbation is adopted for the following analyses.

Difference between R− and Rg-binned snails
As illustrated in Li (2021), the vertical phase snails binned by the Galactocentric radius (R) and the guiding center radius (R g ) are different.The R g -binned snails show higher clarity than the R-binned snails, showing more snail wraps, especially in the outer region of the phase space.This is because the R g -binned stars have similar average orbital radii, and thus they feel similar vertical potential that leads to the winding of the phase snail.On the other hand, in the Rbinned snail, stars have a wide guiding center radius distribution.The R-binned snail can be thought as a combination of different R g -binned snails.Therefore, this R g -mixing effect blurs the R-binned snail.This difference is clearly seen in Fig. 2, in which Rand R g -binned snails at two different radii from the test particle simulation are compared.Li (2021) found that the quantitatively measured shapes of Rand R g -binned snails are similar.However, in order to derive the vertical potential with the snail intersections, we need to carefully quantify their difference.As an intuitive thought, stars of both snails should feel the same potential.In Fig. 2, we overplot the toy model snails as the blue solid lines, which are built by the 1D vertical orbit integration under the vertical potential at Ψ(z)| R=8 kpc (left two panels) and Ψ(z)| R=9 kpc (right two panels), ignoring the radial and azimuthal motions.The model snails initially are energy grids and located at the negative V z axis, mimicking the impulse approximation with a negative velocity kick.For a star with a given vertical energy, we can transform it between z − V z space and Ω z − θ z space, or conversely.We evolve the initial energy grids 500 Myr as the test particle simulation with the 1D vertical orbit integration, which then naturally wrap into the toy model snails.As shown in Fig. 2, especially in the third panel, these toy model snails seem to follow the ridges of the R-binned snails well, especially in the inner region with lower vertical energies.On the other hand, in the very outer regions of the phase space, the toy model snails deviate from the snail density maps in the simulation, which is due to the ignoring of the in-plane motion.Nevertheless, this toy model snail derived from the 1D phase-space orbit integration works well in the inner region, where we can also clearly detect the phase snail in observations.However, on the R g -binned snail density maps, we find that these toy model snails, i.e. the blue lines in Fig. 2, lie in the outer edges of the snails.It is more clearly illustrated in the inner radius of R g = 8 kpc, i.e. the second panel in Fig. 2.This confirms that the Rand R g -binned snails are quantitatively different.

Radius correction for the Rg-based phase snail
Figure 2. The difference between Rand Rg-binned snails.The left two panels are for R = 8 and Rg = 8 kpc, while the right two panels are for R = 9 and Rg = 9 kpc, as labelled on the top left corners.The bin width is 0.6 kpc.All the maps are colored by the number density.The model snail shape curves from the phase-space orbit integration of the vertical potential grids (initially located at Vz < 0 axis) are plotted by the solid lines.The blue lines are for the R-binned snails at the same radii.The red lines are from the potential radii corrected for the Rg-binned snails, with potential radii corrected from Fig. 3, which is about 0.14 kpc larger than the bin centers.Notice that, the blue lines are usually located in the outer edges of the Rg-binned snails, while the red lines are located close to the ridges of the snails.This verifies the necessity of potential radius correction for a Rg-binned snail.
Considering a pure circular orbit with a zero radial action J R = 0, R g is exactly the Galactocentric radius R. For a star with a non-zero radial action, its orbit is rosette-like with radial oscillations around the guiding center radius.The star spends more time in the outer region of the guiding center radius because of the smaller velocity near the apocenter.Therefore, on average, the star feels a vertical potential at a radius slightly larger than R g .We call this time weighted average orbital radius as the vertical potential radius R pot .
In principle, this potential radius can be estimated under the Milky Way's potential for a star with a given radial action, which also varies for different stars.Thus, the calculation is model dependent and time consuming.Here, we propose an empirical approach to derive this potential radius based on the distribution of the guiding center radius of a stellar sample.Considering a narrow Galactocentric radial bin centered at R, the stars inside this bin should on average feel the vertical potential at R. The guiding center radius distribution of these stars would be skewed to a smaller radius.This is similar to the asymmetric drift effect.This Galactocentric radius R can be considered as the potential radius R pot for these stars.Taking the medians of the R g distributions of stars at different radii, we can build an empirical relation between R g and R pot .We then utilize an interpolation of this relation to relate a R g -binned snail to a certain potential radius R pot .
In the test particle simulation, for particles located in different R bins, the difference between the median of R g distribution and the Galactocentric radial bin center (i.e.R pot ) is shown in Fig. 3.The bin width is 0.6 kpc, which is the same in the following simulation analyses.Each Galactocentric radial bin usually corresponds to a R g distribution covering about 2-3 kpc.As a result, the median R g is about 0.1-0.15kpc smaller than the bin center radius, i.e. the potential radius R pot .
To verify this empirical correction, we overplot the toy model snails at R pot on Fig. 2 as the red lines.As shown in the second and fourth panels, the red lines are located closer to the ridge of the snail than the blue lines.They also stay in the inner regions of the blue lines, indicating a smaller or shallower vertical potential.Note that, in the very outer region, the model snails also deviate from the snail density maps, indicating that 1D orbital integration is not proper at such vertical height (> 1.5 kpc).
Ignoring this correction will result in a systematic underestimation in measuring the vertical potential from the R gbinned phase snail, since the vertical potential profile at R pot is slightly shallower than that at R g .This effect could be non-negligible at smaller radii, where the potential changes significantly with a small radius difference.In addition, this test particle simulation imitates a thin disc tracer population.For the observations, the difference between the potential radius R pot and R g could be larger.The correction also depends on the sample selection, such as selecting stars with colder orbits will lead to a smaller difference between R pot and R g .

Refined linear interpolation method for snail intersections
As illustrated in Fig. 1, our method relies on the correct interpolations of the Z max and V z,max values for the intersections between the phase snail and the V z axis and z axis, respectively.A straightforward approach for the interpolation is to directly use the Z max or V z,max values of the adjacent intersections.For example, Z max of point 4 can be estimated with the average Z max of points 3 and 5, and V z,max of point 3 can be estimated with linear interpolation of V z,max of points 2 and 4.However, in order to derive the vertical potential with less systematical errors, physical quantities linearly related to the vertical frequency Ω z are better choices for the interpolation.We test different parameters using the model snails with the results shown in Fig. 4. The model snails at R pot = 8 kpc and R pot = 8.14 kpc (corresponding to the snail at R g = 8 kpc) are displayed as the blue and red solid lines in the left panel of Fig. 4, respectively.The intersections are defined by the action angles where θ z = k 2 π (k is an integer).They are marked as dots in the left panel, and their corresponding Ω z are indicated by the vertical dashed lines in the middle panels of Fig.

4.
For the z-axis intersections, we could utilize V z,max or E z = 1 2 V 2 z,max as the linear interpolation quantity.The relations of Ω z with E z and V z,max are shown as the solid and dashed lines in the second panel of Fig. 4, respectively.Apparently, E z has a non-linear dependence on Ω z , while V z,max has a roughly good linear relation with Ω z , in the region considered.Therefore, using E z as the linear interpolation quantity will result in an overestimation, while utilizing V z,max can significantly reduce the systematic error in the interpolation.For this reason, we adopt V z,max as the linear interpolation quantity to derive Ψ(Z max ) = E z for the z-axis intersections.
For the V z -axis intersections, we explore the relation between Z max and Ω z , shown as the solid lines in the third panel of Fig. 4. It is also a curved relation, similar to the relation between E z and Ω z .Thus, using Z max directly as the linear interpolation quantity will result in an overestimation.As there is squared relation between E z and V z,max , we explore the relation between the square root of Z max , i.e. √ Z max , and Ω z , illustrated by the dashed lines in the third panel of Fig.

The relation between
√ Z max and Ω z seems linear in the region covered by the vertical dashed lines.Thus, √ Z max is adopt as the linear interpolation quantity to derive Z max for the V z -axis intersections.
Though the above two approximate linear relations are derived under the potential of Irrgang et al. (2013), they should remain valid in real observations for two reasons.First, the shapes of popular Milky Way potentials are similar after shifting to the same Solar radius.Second, in other radial bins, there are also good linear relations of Ω z with √ Z max and V z,max .The linear interpolation of these two quantities should work well at different radii.As a summary, we re-fine our interpolation method for the (i + 1)-th intersections as following: (6) Z max,i and V z,max,i are the absolute values measured for the intersections with phase space coordinates of (z i , V z,i ) on the z-axis and V z -axis, respectively.
Combining with Eq. 3, we can obtain the (Z max , Ψ(Z max )) for N-2 intersections (excluding the first and last snail intersections), where N is the number of intersections detected for a snail.The interpolation is performed for the intersections measured from the model snails, and the resultant potential data points are shown as the filled symbols in the right panel of Fig. 4.They have a good consistency with the input potential profiles, with residuals quite close to zero, as illustrated by the open symbols.However, ignoring the difference between R g and R pot will result in a systematic underestimation of up to 5%, as indicated by the black open symbols in the top right panel of Fig. 4.
According to the 1D model snail test, the interpolation method defined by Eq. 6 is accurate, showing much less systematic uncertainty.In real observation, the main error source in the observation should arise from the determination of these intersection positions, and the ignorance of the difference between the vertical potential profiles at R g and R pot .The derived (Z max , Ψ(Z max )) potential points provide a direct and model independent measurement of the vertical potential profile of our Milky Way.It can be compared with the popular Milky Way potentials, and applied with a further mass modelling to construct the Galactic mass distribution.

Intersections measurement and uncertainties
This section describes our method to measure the snail intersections with both axes.We start with the test particle simulation.Although the R-binned snails on average feel the vertical potential requiring no potential radius correction, they are more dispersed than the R g -binned snails.In order to reduce the errors of intersection measurements, we prefer to measure the intersections based on the R g -binned snails and apply the interpolation method, which will also be performed on the observational data in the next section.
We take five radial bins centered at R g = 7, 8, 9, 10, 11 kpc with the bin width of 0.6 kpc.The phase space density distributions of the snails are shown in the left column of Fig. 5, with a pixel size of 0.026 kpc × 1.8 km s −1 .In observations, the snail number density map is usually convolved with a Gaussian kernel to get a smooth background (Laporte et al. 2019;Li & Shen 2020).The density contrast map ∆N after subtracting the smooth background significantly enhances the snail signal.We also show the density contrast maps (with a Gaussian kernel dispersion of σ = 3 pixels) for our R g -binned snails in the second column of Fig. 5.
In order to measure the intersections, we extract the number density profiles along two stripes at V z = 0 (along z axis) and z = 0 (along V z axis), with widths of 20 km s −1 and 0.3 kpc, respectively.These widths are similar to the widths of the snails in the z and V z directions.Similar to the 2D density contrast map, we perform 1D Gaussian kernel convolution to the extracted number density profiles.The density contrast profiles are displayed in the third and fourth columns of Fig. 5, for the two stripes at V z = 0 and z = 0, respectively.In these 1D number density contrast profiles, there are several peaks related to the locations of intersections.By visual inspection, we select the rough regions where peaks emerge.We ignore several peaks located in the inner region of the phase space, where the peaks could be blurred by the high density background in observations.For these regions, we fit the density contrast profiles with a Gaussian function with the peak corresponding to the location of an intersection.The best-fit Gaussian profiles are shown as the red lines in the right two columns of Fig. 5, while the intersections are indicated by the light green dashed lines.The error bars of the intersections are taken from the maximum of the error from the Gaussian fitting and the pixel size, of which the former is usually smaller than the latter.The measured intersections are also marked as the light green dots in the original and contrast density maps, shown in the left two columns of Fig.

5.
For the measured intersections, we can apply the interpolation method in Eq. 6 to obtain vertical potential measurements at several vertical heights.The results are shown in Fig. 6.In the lower panels, the interpolated points are well consistent with the vertical potentials at the corresponding potential radii, i.e. the red lines, which are derived from the relation displayed in Fig. 3.As explained previously, these points show systematic underestimation, compared to the potential just at the guiding center radii (the blue lines).This underestimation is more significant at smaller radii.The residuals between the derived potential data points and the potential profiles at R pot and R g are indicated by the red and blue points in the top panels of Fig. 6, respectively.Compared to the potential at R pot , the interpolated potential points are consistent with deviations less than 5%.However, compared to the vertical potential profiles at R g , they manifest systematic underestimation with slightly shallower profiles, up to ∼ 5% at the inner most radius.Note that, our test particle simulation mimics a cold thin disc, which shows only about 0.1-0.15kpc difference between R g and R pot .In the observations, this potential radius correction could be larger, resulting in a more significant underestimation.As shown by the error bars in the upper panels of Fig. 6, the uncertainty due to the intersections measurement is about 10%.As we do not include background particles in the measurement, this uncertainty could also be larger in the observations.

DATA
We try to obtain the disk vertical potential measurements with the interpolation method of the phase snail intersections based on Gaia Data Release 3 (Gaia Collaboration et  al. 2023).Gaia DR3 RVS provides radial velocity measurements for about 33.8 million objects.We take the geometric distances from Bailer-Jones et al. ( 2021), which are based on a Bayesian approach that uses a prior constructed from a three-dimensional model of the Milky Way.We apply a distance quality cut, requiring the relative parallax error (σ ϖ ) smaller than 20 per cent (i.e.ϖ/σ ϖ > 5).In terms of the astrometric measurements, we require the renormalised unit weight error (RUWE) to be smaller than 1.4.The final sample contains 27,209,191 stars.
The coordinates system adopted in this paper is the classical Galactocentric cylindrical system (R, ϕ, z), with the Galactocentric radius R increasing radially outward, ϕ increasing in the direction of Galactic rotation, and positive z pointing toward the North Galactic Pole.The Sun is located at (−8.34, 0, 0.027) kpc in the Galactic Cartesian coordinates system (Chen et al. 1999;Reid et al. 2014).The Solar motion relative to the local standard of rest (LSR) is (U ⊙ , V ⊙ , W ⊙ ) = (7.01,12.09, 7.01) km s −1 (Huang et al. 2016), and the circular velocity of LSR is adopted as 240 km s −1 (Reid et al. 2014).

Measuring the snail intersections in the observation
Although R-binned snails directly represent the vertical potential, they are more dispersed than R g -binned snails.In order to reduce the uncertainty in measuring the snail intersections, we utilize R g -binned snails in observations.The guiding center radius R g is calculated as the radius of a circular orbit with the same angular momentum (L z ) as the star, i.e.R g = RV ϕ V c (R g ), where V ϕ is the azimuthal velocity and V c is the rotation curve.We utilize the rotation curve from Huang et al. (2016), with a circular velocity of 240 ± 6 km s −1 at the Solar radius (R 0 =8.34 kpc), to build a relation of L z − R g .By mapping the observed angular momentum of a star to the derived L z − R g relation, we can obtain its R g through an interpolation.We explore nine radial bins with centers at R g = 7.5, 8., 8.34, 8.5, 9, 9.5, 10, 10.5, 11 kpc, with a bin width of ∆R g = 0.5 kpc.More inner radii are not explored in this work as the snails could be more complicated, while at the outer radii, measurements could be dominated by the noise.We trace the snails in a phase space region of |z|< 1.5 kpc and |V z | < 75 km s −1 , as the outer region is dominated by the numerical noise.
The phase space density distributions of R g -binned snails are shown in the left column of Fig. 7, with their density contrast maps displayed in the second column.Note that, compared to the simulated snails, the observed ones contain large fractions of background stars.Nevertheless, we have one order of magnitude larger sample size.The left two columns show clear radial variation of the phase snails, as snails at outer radii are more elongated along z axis.This is due to the fact that the outer radius has shallower vertical potential, which results in a larger Z max for a given vertical energy.The density contrast profiles along z and V z axes are shown in the right two columns of Fig. 7.The snail intersections are measured in the same way as in the test particle simulation, and are illustrated by the vertical dashed lines in the right two columns and the light green points in the second column of Fig. 7.
The empirical relation between the guiding center radius and the potential radius for the observational sample is shown in Fig. 8.The corresponding angular momenta (L z ) and the potential radii (R pot ) for the R g bins used are listed in Table 1.The bin width of each R bin is 0.5 kpc, which is the same as the R g bin width.In the region of 6 < R < 10 kpc, the guiding center radius is about 0.4 − 0.5 kpc smaller than the potential radius.This difference is significantly larger than that in the test particle simulation, which could thus result in a significant underestimation if ignoring this R pot correction.At R ∼ 11 kpc, the difference reaches a minimum of about 0.2 kpc.In the outer (R > 12 kpc) region, the relation shows a rapid decreasing pattern, which could be due to insufficient stars.The relation in the inner region is not shown due to the truncation of the rotation curve of Huang et al. ( 2016) at about 5.5 kpc.Note that, this R g vs. R pot relation is empirical, and should change with the sample selection.For example, if we choose stars in colder orbits, the difference between R g and R pot would be smaller.However, such choice would reduce the number statistics to allow for the accurate snail shape measurement.
Based on the measured snail intersections from the observation, we apply the interpolation method as defined by Eq. 6 to obtain Z max and Ψ(Z max ) = 1 2 V 2 z,max of these intersections.The derived snail potential data points are shown in Fig. 9 as the red dots.These data points provide a model independent assessment of Milky Way's vertical potential, and can also be applied to the further mass modelling.

Comparison with the popular Milky Way potentials
We compare our measurements of the Milky Way potential with three popular models of it, which are the Model I from Irrgang et al. (2013) (Irrgang13I), MWPotential2014 from Bovy (2015), and McMillan17 from McMillan (2017), shown as the blue, red and green lines, respectively, in Fig. 9.As these potentials have different Solar positions, they are radially shifted to the same Solar radius, i.e.R 0 = 8.34 kpc.Besides the vertical potentials at the centers of R g bins (solid lines), the potentials at the corresponding potential radii are also displayed with the dashed lines in Fig. 9.The potential radii are calculated from the empirical relation shown in Fig. 8.
The potential points derived from the snail intersections agree very well with the other popular Milky Way potentials    from 7.5 to 11 kpc at |z| < 1 kpc.The interpolated snail intersections are slightly shallower than the vertical potential profiles at R g , and are more consistent with those profiles at the corresponding potential radii.This phenomenon is more obvious at inner radii, where the vertical potential difference between R g and R pot is larger.At outer radii, as the radial gradient of the potential becomes shallower, the vertical potential profiles at R g and R pot are similar to each other.For the current snail intersection measurements, we only explore low vertical regions, which are only |z| ≲ 1 kpc as shown in the second column of Fig. 7.In addition, we have no interpolation for the outermost intersections in the phase space, which would result in smaller regions of the interpolated snail potential points.
Note that, at the inner two radii, i.e.R g = 7.5, 8.0 kpc, there are two data points at z ∼ 0.7 kpc showing significant deviations from the Milky Way potential profiles.This is likely due to the Z max estimation uncertainties for the intersections at z ∼ 0.7 kpc.As shown in the second panel of the top row in Fig. 7, the snail density contrast map at z ∼ 0.7 kpc seems to show two adjacent wraps, with the outer wrap deviating from the snail pattern of the inner wrap.Antoja et al. (2023) also found an additional wrap at V z ∼ 50 km s −1 in the volume at R = 8 kpc, while our additional snail wrap is more clearly seen at z ∼ 0.7 kpc and V z ∼ 0 km s −1 .The outer wrap is more obvious, which results in a larger Z max for that snail intersection.The interpolation transfers this overestimation to the two adjacent snail intersections on V z axis.If we take this Z max overestimation into consideration, the relevant inter-polated snail intersections will become closer to the dashed lines in Fig. 9.

Modelling the global disc
Based on the direct measurements of the vertical potential from the interpolated snail intersections, we perform mass modelling to these potential data points, to constrain the mass distributions of different Galactic components.We apply a mass model consisting of double exponential discs and an NFW (Navarro et al. 1996) dark matter halo.For simplicity, we regard the two discs as very flattened systems, and integrate the density profile along the z direction twice to directly obtain their reduced vertical potential.See the detailed derivation in Appendix A. The total reduced vertical potential at (R, z) is: ) where r = √ R 2 + z 2 and G is the gravitational constant.M d1 , R d1 and z d1 are the total mass, the scale length and scale height of the thin disc, respectively, while M d2 , R d2 and z d2 are for the thick disc.M h and a h are the scale mass and scale radius for the dark matter halo.We use the Markov Chain Monte Carlo (MCMC) technique to fit the interpolated snail intersections and obtain the parameters' estimation.The MCMC package we use is EMCEE (Foreman-Mackey et al. 2013).As our snail potential data points are at low vertical regions, we fix the scale heights of discs as z d1 = 0.3 kpc and z d2 = 0.9 kpc.We also apply Gaussian priors on the radial scale lengths, which are 2.6 ± 0.52 kpc and 3.6 ± 0.72 kpc for the thin and thick discs, respectively (Bland-Hawthorn & Gerhard 2016;McMillan 2017).
In the mass modeling, we choose the corrected potential radius for each phase snail to derive the vertical potential Ψ(R pot , z). Results of the global mass modeling of the interpolated snail intersections are shown in Fig. 10.For comparison, the potential profiles from MWPotential2014, and the profiles at the uncorrected radius (i.e. the center of each R g bin) Ψ(R g , z) from our best-fit potential model are also displayed.The model predicted vertical potential at R pot , i.e. the black dashed lines, is well consistent with the data points.The vertical potential profiles at inner radii are consistent with those from the MWPotential2014 (red lines), for both R g and R pot .At outer radii (R g > 10 kpc), the snail intersections prefer steeper vertical potential, resulting in an overestimation of the outer vertical potential.Although we As mentioned in Section 2, the snail at R g = 8.34 kpc can not represent the vertical potential at the Solar neighborhood.According to the empirical relation shown in Fig. 8, a potential radius of R pot = R 0 = 8.34 kpc is related to the snail at R g = 7.91 kpc, showing a 0.43 kpc difference.Besides, according to the top three rows in Fig. 7, snails seem to show one more wrap in the outer phase space.Therefore, in order to study the vertical potential from the snail intersections in the Solar neighborhood, we explore the snail within |R g − 7.91| < 0.25 kpc and measure the intersections for one more wrap than that in the previous Section.The measurement of the snail intersections is shown in Fig. 11, which is similar to Fig. 7.There are 12 snail intersections, resulting in 10 snail potential data points up to z ∼ 1.2 kpc after the interpolation.
The interpolated snail intersections are shown in the upper panel of Fig. 12 as the red dots.For comparison, we also show the Solar vertical potential from three popular Milky Way potentials (i.e.Irrgang13I, MWPotential2014, McMil-lan17) as the colored lines.The snail intersections are well consistent with the three Milky Way potentials, except for two points at z ∼ 0.7 kpc.These two points could suffer from the same problem as explained in Section 4.2 for the snail at R g = 7.5 kpc.
We apply a simple mass model to the snail intersections at R pot = R 0 , which consists of an exponential thin disc, a razor thin gaseous disc and a constant local dark matter density (ρ dm ).This model is the same as that in Guo et al. (2020).The vertical potential is therefore written as: (8) where Σ ⋆ and Z h are the total surface density and the scale height of the stellar disc, respectively.The total gas surface density Σ gas is set as 13.2 M ⊙ pc −2 (Flynn et al. 2006).The free parameters left are p = (Σ ⋆ , Z h , ρ dm ).See the derivation of this potential model in Appendix B or in Guo et al. (2020).
To consider the uncertainty in the Z max estimation, we perform 1000 Monte Carlo realizations according to the error of each intersection.The uncertainty of Ψ(Z max ) is absorbed into the following log-likelihood: where Ψ obs,i and σ Ψ,i are the interpolated vertical potential and its uncertainty for the i-th snail intersection, while Ψ model ,i (Z max,i ) is the model vertical potential at the related Z max,i for this intersection.We further apply Gaussian priors for Σ ⋆ and Z h , which are 37.0 ± 5.3 M ⊙ pc −2 (Guo et al. 2020) and 0.5 ± 0.1 kpc, respectively.For the local dark matter density, we just constrain it in a range of [0, 0.05] M ⊙ pc −3 , with a flat prior.
The posterior distributions of these parameters from the MCMC method are shown in Fig. 13.The medians and 1σ uncertainties of the 1D marginalized distributions are 33.1 ± 5.2 M ⊙ pc −2 , 0.612 ± 0.085 kpc, 0.0150 ± 0.0031 M ⊙ pc −3 for the disc surface density, the disc scale height and the local dark matter density, respectively.The best-fit parameters are in coincidence with the medians, as shown by the dashed lines in the histogram panels in Fig. 13.
The model predicted vertical potential profile and the corresponding total surface density profile are shown as the black solid lines in Fig. 12, with the shaded regions representing their 1σ errors.The model predicted vertical potential is still systematically smaller than the three Milky Way potentials.This is mainly due to the large uncertainties of the two snail intersections at z ∼ 0.7 kpc (resulting from the possible additional wrap), rather than due to the potential radius correction.For the total surface density, i.e. the slope of the vertical potential, the model predicted one shows significant difference compared to the three Milky Way potentials.There are two reasons.The first is that we set the constant gas surface density.This would influence the surface density profile at the low vertical region (z < 0.2 kpc), as the scale height of the neutral hydrogen is about 0.085 kpc (McMillan 2017).The second reason is the possible overestimation of Z max for the two snail intersections at z ∼ 0.7 kpc, which will decrease the potential slope at the lower vertical height and increase the potential slope at the higher vertical height.The shape of the total surface density profile is then bent downward accordingly, which may explain the difference at 0.5 < z < 1.0 kpc.

Comparison with previous works
There are a lot of works trying to measure the vertical potential and mass distribution in the Solar neighborhood (e.g.Zhang et al. 2013;Xia et al. 2016;Hagen & Helmi 2018;Schutz et al. 2018;Sivertsson et al. 2018;Buch et al. 2019;Guo et al. 2020;Salomon et al. 2020).These works utilize dynamical methods, such as the distribution function and the vertical Jeans equation, to model the stellar kinematics, based on the equilibrium state assumption.Here, we compare the vertical potential profile and the total surface density profile of our independent measurement with those from previous works, e.g.Holmberg & Flynn (2000); Korchagin et al. (2003); Holmberg & Flynn (2004); Bovy & Rix As shown in the upper panel of Fig. 14, our model predicted vertical potential profile, i.e. the median of all potential profiles derived from the posterior parameter sets, shows quite good agreement with the profiles from Li & Widrow (2021) and Holmberg & Flynn (2000).However, consider-ing the possible overestimation of Z max at z ∼ 0.7 kpc, the interpolated snail intersections would be more similar to our previous work (Guo et al. 2020).The profiles from Zhang et al. (2013) seem to slightly overestimate the vertical potential at z ≲ 1 kpc, relative to other works.In the lower panel, if ignoring the constant gas surface density applied, our surface density profile is best consistent with that from Li & Widrow (2021).However, as explained previously, the shape in the region of about 0.5 < z < 1.0 kpc should be bent upward slightly.Thus the intrinsic Σ tot profile might be more similar  to our previous work (Guo et al. 2020).The profiles from Zhang et al. (2013) seem to be slightly larger in the lower region (z < 1 kpc) and smaller in the higher region (z > 1 kpc), where Piffl et al. (2014) also shows a shallower profile.Taking the constant Σ gas and the probable overestimation of Z max at z ∼ 0.7 kpc into consideration, our model predicted surface density profile could be consistent with all the measurements shown by the different symbols in the lower panel of Fig. 14.
For the local dark matter density, our model predicts ρ dm = 0.0150±0.0031M ⊙ pc −3 , which is consistent with our previous measurement of 0.0133 +0.0024 −0.0022 M ⊙ pc −3 (Guo et al. 2020) and other works within the error bars.It is larger than some other works (Read 2014;de Salas & Widmark 2021).However, if we consider the probable overestimation of Z max at z ∼ 0.7 kpc, the Σ tot profile should be bent upward in the region of 0.5 < z < 1.0 kpc.It means that the total surface density will increase more quickly in the lower region, and more slowly in the higher region.Therefore, the scale height of the stellar disc should be smaller, which is also more reasonable as the model predicted Z h = 0.612 ± 0.085 kpc is larger for the thin disc (McKee et al. 2015).The local dark matter density should be slightly smaller considering the smaller slope of Σ tot in the higher region, where the constant dark matter dominates the profile.Comparing the model predicted vertical potential (top panel) and the total surface density profile (bottom panel) at Rpot = R0 with previous results in the literature (as labelled).The black lines show our modeling results by fitting the interpolated snail intersections (red points), with shaded regions indicating the 1σ errors.Our results best agree with the profiles from Holmberg & Flynn (2000) (black dashed lines) and Li & Widrow (2021) (blue solid lines).

Comparison with the phase snail modeling in previous works
Since the discovery of the phase snail, many works have tried to measure the Milky Way mass distribution and potential with the phase snail.In the series works of Widmark et al. (2021aWidmark et al. ( ,b, 2022a,b),b), the authors constructed the analytical snail model through the 1D phase-space orbital integration and a sinusoidal function of the action angle.Assuming a constant amplitude relative to the background, they obtained an estimation of the vertical potential embedded in the phasespace orbital integration by fitting the 2D phase space distribution in the observation.In their model, the orbital integration requires a potential model and the time since the perturbation.In two other works, Li & Widrow (2021, 2023) simultaneously modeled the gravitational potential and phase space distribution function of giant stars under the steadystate assumption of the disk, with clear snail features shown in their residual maps.In our previous work (Guo et al. 2022), we applied a geometric snail model superimposed on a smooth background in the phase space to fit the north-south asymmetric vertical velocity dispersion profile.The phase space distribution of the smooth background component is built through the vertical Jeans equation and the background potential, while the snail model is constructed empirically from both the observation and test particle simulations.In a different work, Price-Whelan et al. ( 2021) utilized Orbital Torus Imaging, which makes use of kinematic measurements and element abundances, to investigate the metallicity distribution in the phase space.Under the steady-state, stellar labels (e.g.metallicity) should vary systematically with orbit characteristics such as actions, while keeping invariant with respect to the action angle (θ z ).By minimizing the residual metal distribution as a function of θ z , they also obtained an estimation of the Milky Way mass model.
There are many works on the phase snail focusing on the measurement of the time since the perturbation, which is found to vary with the angular momentum (Antoja et al. 2023;Darragh-Ford et al. 2023).Antoja et al. (2023) measured the snail shapes first, which vary with the angular momentum and the azimuth significantly.They then inferred the time since the perturbation by the difference of vertical frequencies of snail intersections.The derived impact time varies with both the angular momentum and the vertical height.In their work, a Milky Way potential model is needed to derive the vertical frequencies of snail intersections, though different potential models lead to similar results.In Darragh-Ford et al. (2023), the time since the perturbation is measured from the slope of the action angle (θ z ) as a function of the vertical frequency (Ω z ).They also need a potential model to derive θ z and Ω z of each star, which are also utilized to extract the snail shape.In order to find a suitable potential, they ran the action calculations in a grid of potential parameter values, and visually picked a potential model in which the snail is close to a straight line in the space of √ J z − θ z .In our work, the toy model snails in Section 2.2 are built through the 1D orbital integration similar to Widmark et al. (2021a).However, these toy model snails are only used to distinguish the difference between Rand R g -binned snails and to refine the interpolation method.Different from the previous works, the linear interpolation method here does not need to assume any potential model.In addition, although the guiding radius calculation is based on a rotation curve, which is an implication of the underlying mass distribution, we can actually separate the sample according to the angular momentum, which is directly derived from the observation.In a summary, our method provides a new way to directly obtain potential measurements, independent from the mass modeling and the assumption of the time since the perturbation.

Connection between the phase snail and the
North/South asymmetry In our previous work (Guo et al. 2022), we found that the north and south asymmetry of the Galactic plane can be quantitatively explained by a simple empirical snail model.The number density difference profile between north and south, i.e.A z = (ν(z) − ν(−z))/(ν(z) + ν(−z)), shows a wavelike fluctuation (Widrow et al. 2012;Bennett & Bovy 2019), which is caused by the snail intersecting with the z axis at different vertical heights.To inspect this effect and check if it could be utilized to measure the snail intersections, we investigate the number density difference profiles between the positive and negative vertical height (A z ) and vertical veloc- We select the same stripes along z and V z axes as Fig. 7 to check the number density difference profiles, in order to avoid the cumulative effect of different snail wraps when using the whole phase space.Note that, we ignore the selection effects and regard they are the same at positive and negative z or V z .The results are shown as the black dots in Fig. 15.The snail intersections measured in the previous section are indicated by the vertical dashed lines.We also show Gaussian smoothed profiles for the data points, shown as the black dashed lines in Fig. 15.
For most snail intersections, they are consistent with the locations of peaks and troughs, especially at smaller vertical heights and velocities.This is due to the higher number density for regions closer to the phase space origin.At larger vertical heights and velocities, due to the large numerical noise, the peaks and troughs are sometimes not obvious.It results in a weak correlation with the snail intersections.Besides, the snail wraps are tighter at inner radii, which results in smaller absolute difference between the adjacent intersections along z and V z axes.Therefore, the adjacent peaks and troughs of A z and A Vz profiles could be mixed up, as indicated by the outer two intersections in the top left panel of Fig. 15.Another interesting result is the amplitudes of the A z and A Vz profiles, which vary with the guiding center radii, the vertical heights and vertical velocities.The amplitudes may include the information about the perturbation strength, which will be investigated in our future work.One important assumption in our work is that the vertical motion is separable from the in-plane motions and the vertical energy is conserved, which thus relates Z max and V z,max of snail intersections as Eq. 3.This is an appropriate assumption at lower vertical heights, but breaks down at high vertical regions.As shown in Fig. 2, our toy model snail from 1D vertical integration can roughly trace the ridge line of the snail in the test particle simulation.At the outer region of phase space (about z > 1.5 kpc), the snails always stay inside the toy model snails, due to the ignorance of the in-plane motions.To trace the vertical potential from snails at large vertical heights, 3D orbit integration is needed.Otherwise, it would result in a shallower profile at high regions.

Some caveats in our method
Another caveat is the empirical potential radius correction for R g -binned snails.Though our empirical correction method is verified with the test particle simulation, it should be carefully applied in the observations.Both the sample selection and the data binning may have slight influence on the resultant correction relation.The decline at the tail in Fig. 8 could be due to the incompleteness at large radii.Besides, stars in cold orbits and hot orbits should have different correction relations, as their radial actions are different.As illustrated in Li & Shen (2020), the hotter orbits phase-mix more quickly and may show no snail features in the phase space.We select all stars in this work, which results in a 0.4 − 0.5 kpc difference between R g and R pot .If only stars with smaller radial actions are selected, this difference could be smaller.The accurate correction is difficult to perform for each star, which may need the 3D orbit integration.More detailed studies are needed to accurately quantify the connection between the shapes of R g snails and the underlying vertical potential.

The self-gravity effect
The self-gravity is ignored in our work, especially when testing our method with the test particle simulation and building the model snails.However, the self-gravity is important in modelling the evolution of the snails, especially the tightness of the snails.The phase snails found in N-body simulations are not as tightly wound as in the observations (Laporte et al. 2019;Bennett et al. 2022).Darling & Widrow (2019) compared N-body simulations of a test particle disc with fully self-consistent ones.They found that the snails in the simulations with the self-gravity are less tightly wound and more difficult to discern.More recently, Widrow (2023) explored the influence of the self-gravity on the amplitude and pitch angle of the phase snails within the framework of the shearing box approximation.With the self-gravity, the surface density is amplified due to the swing amplification, which thus sets off new vertical phase snails.The snails are then less tightly wound while stronger than those in the case without the self-gravity.The slope in the frequency and ac-tion angle space therefore indicates a smaller perturbation time, related to the time between the swing amplification peak and the snapshot time.
Our interpolation method bases on the linear relations between quantities and the vertical frequency.When the self-gravity is considered, the vertical frequencies of snail intersections are gradually changed, with different extent.Though the phase angle difference between adjacent intersections is still π/2 and thus the vertical frequency difference keeps constant, the approximate linear relations between quantities and Ω z might be changed.On the other hand, these relations are generally accurate for different vertical potential profiles at different radii.As a balance, it is difficult to assess if the interpolation of averaging the adjacent intersections would be problematic or not.In addition, the self-gravity would also gradually change the vertical density profile and thus the vertical potential profile.The shape of the phase snail contains both the information of the vertical potential profile before the perturbation and its gradual variation due to the self-gravity.The vertical potential traced by the snail seems neither to be the profile before the perturbation, nor the current profile.Finally, the strength of the perturbation, i.e. the fraction of stars being strongly perturbed or in the snail structure, is important in assessing the influence of the self-gravity.Studies on the amplitude of the phase snail (e.g.Guo et al. 2022;Alinder et al. 2023;Frankel et al. 2023) are helpful to understand it.All these considerations need more detailed studies with the self-gravity considered.We leave them to the future work.

CONCLUSION
The vertical phase snail is a direct sign of dis-equilibrium of our Milky Way, which is the result of incomplete phasemixing.Nevertheless, the snail wraps mainly under the vertical potential, which means it contains the information of the Galactic mass distribution of different components.In this work, we propose a novel method to measure the Galactic vertical potential, utilizing the intersections between the phase snails and z/V z axes (i.e. the turning-around and midplane points).The vertical motion is assumed to be decoupled with the in-plane motion, which is a proper assumption at low vertical height (z < 1.5 kpc).The snail intersections have known maximum vertical heights (Z max ) for the turningaround points or vertical velocities (V z,max ) for the mid-plane ones.Applying an interpolation method, we can thus obtain (Z max , 1 2 V 2 z,max ) for each intersection, which then provides a direct and model independent measurement of the vertical potential.
This method has been tested utilizing a test particle simulation.We find approximate linear relations of V z,max and √ Z max with the vertical frequency Ω z .Thus, they can be utilized in the linear interpolation between the adjacent snail intersections, which have a constant action angle interval θz = π/2 and thus a constant ∆Ω z .In addition, we find that the snail binned by the guiding center radius is quantitatively different from that binned by the Galactocentric radius.Direct shape comparison shows that the R g -binned snail stays slightly inside the R-binned snail, which means it stands for a vertical potential shallower than that for the R-binned snail.This is due to the fact that stars on a non-circular orbit stay longer outside the guiding center radius.Thus on average, the R g -binned snail stands for a vertical potential with a potential radius larger than R g , i.e.R pot > R g .This difference can also be understood as the asymmetric drift: the median R g of a sample in a given radial range is smaller than the center of the radial bin, i.e.R pot .The difference between R pot and R g depends on the radius, the radial action of the tracer and thus the sample selection.We apply an empirical method to relate R g with R pot .From our test particle simulation, the difference between R pot and R g is about 0.1-0.15kpc, which results in about 5% systematic underestimation of the vertical potential with R g not corrected.After the correction from R g to R pot , our interpolation method obtains an estimation of the vertical potential with less than 5% random errors.
We apply our method to the Gaia DR3 for several R g bins.The empirical relation between the potential radius and guiding center radius for our sample shows that R g is about 0.4 kpc smaller than R pot , which needs to be corrected in the potential comparison and the further mass modeling.We apply a Gaussian kernel convolution to obtain the density contrast profiles for stripes along z/V z axes.The snail intersections are measured by Gaussian fitting to visually selected regions where peaks are located.With the measured Z max for the turning-around points and V z,max for the mid-plane points, we apply the linear interpolation method as refined in Eq. 6 to obtain (Z max , 1 2 V 2 z,max ) for these intersections.Note that, at R g = 7.5, 8.0 kpc, the maximum vertical heights of the intersections at z ∼ 0.7 kpc seem to be overestimated, resulting in overestimation of Z max for adjacent snail intersections.In addition, the locations of intersections are also coincident with some peaks and troughs of the difference profiles between the positive and negative vertical heights (A z ) and vertical velocities (A Vz ), which might be utilized to locate the snail intersections.
We compare the potential measurements for the snail intersections with three popular Milky Way potentials, i.e. the Model I from Irrgang et al. (2013), the MWPotential2014 from Bovy (2015) and the McMillan17 from McMillan (2017).Our measurements are more consistent with the vertical potential at the potential radii, which are obviously shallower than the vertical potential at the smaller guiding center radii.At larger radii, the difference is smaller as the radial variation of the vertical potential is smaller.We apply further mass modeling to the snail intersections at different radii, with the scale heights fixed and priors applied to the scale lengths of the double exponential disc.The resultant vertical potential is well fitted to the snail intersections and consistent with the potential from Bovy (2015) at small radii.However, it shows an overestimation at large radii.
From the empirical relation between R g and R pot for our sample, we attempt to measure the vertical potential in the Solar neighborhood, which is related to the snail at R g = 7.91 kpc.We manage to measure the snail intersections for one more wrap with relatively large uncertainties.The resultant potential measurements for the snail intersections are consistent with the Solar vertical potential from the three popular Milky Way potentials.We fit the interpolated vertical potential of these intersections with a mass model consisting of a thin stellar disc, a razor thin gaseous disc and a constant local dark matter density.Under the priors on the stellar total surface density and the disc's scale height, we obtain a local dark matter density of ρ dm = 0.0150±0.0031M ⊙ pc −3 , which is consistent with previous works.The model predicted vertical potential profile and the vertical total surface density profile are consistent with previous results, especially with Holmberg & Flynn (2000) and Li & Widrow (2021).However, if we take the possible overestimation of Z max at z ∼ 0.7 kpc into consideration, the intrinsic vertical potential could be more consistent with Guo et al. (2020).The profile of Σ tot would be bent upwards, resulting in a smaller stellar vertical scale height and a smaller local dark matter density.
In principle, this model independent method can be applied to more radial bins with overlapped stars.In this way, we can better compare with different known Milky Way potentials, and provide a stronger constraint on the Milky Way potential.Our work indicates that the vertical phase snail can indeed provide a new way to study the Milky Way's vertical potential.In our future work, we will try to measure the vertical potential utilizing the whole snail shape curve, rather than just the snail intersections, which might also better constrain the Milky Way potential.

Figure 1 .
Figure 1.Sketch plot of the intersections between the phase snail and z/Vz axes.The blue line is an Archimedean spiral.The horizontal and vertical dashed lines represent the z and Vz axes, respectively.The indices of intersections are also labelled, with smaller indexes for inner intersections.The turning-around (Vz = 0) and mid-plane points (z = 0) are marked in odd and even numbers, respectively.

Figure 3 .
Figure3.The difference between a radial bin center and its median guiding center radius (Rg) for the test particle simulation.The bin width of each radius is 0.6 kpc.This relation is utilized to relate a Rg-binned snail to the potential radius (Rpot ).

Figure 4 .
Figure 4. Vertical potential derived from idealized phase snail calculated from the 1D orbital integration.The blue and red lines represent the snail shape curves at Rpot = 8 kpc and Rpot = 8.14 kpc, corresponding to the snails binned by |R − 8| < 0.3 kpc and |Rg − 8| < 0.3 kpc, respectively.Left panel: snail shape curves in the z − Vz phase space, with snail intersections marked by dots (for z < 1.5 kpc).Four intersections labelled by the larger red squares are marked in all panels to help understand the interpolation method.Middle left panel: the vertical energy (Ez in the left axis, solid lines) and the maximum vertical velocity (Vz,max in the right axis, dashed lines) as a function of the vertical frequency (Ωz) for the two snails shown in the left panel.The vertical dashed lines indicate the frequencies of the intersections.Middle right panel: same as the middle left panel but for the maximum vertical height (Zmax in the left axis, solid lines) and its square root ( √ Zmax in the right axis, dashed lines).Note that within the region covered by the vertical lines, the dashed lines show better linear relations with the vertical frequency than the solid lines.Right panel: snails presented as Ez versus Zmax, i.e. the vertical potential profiles.The filled symbols are derived utilizing a linear interpolation function based on the relations between Ωz and Vz,max and √ Zmax (dashed lines in the middle panels), for the intersections marked in the left panel.In the top panel, the open symbols are the relative residuals between the filled symbols and the model potential curves.The black symbols are the relative residuals between the red filled symbols and the blue solid line, which shows a systematical underestimation for using a Rg binned snail without considering the potential radius correction.

Figure 5 .
Figure 5. Measure the intersections for the simulated phase snails.The left column shows the snails at different guiding center radii (as indicated), while the second column presents their density contrast maps.The pixel size is 0.026 kpc × 1.8 km s −1 .The third and fourth columns display the density contrast profiles (grey lines) for stripes at Vz = 0 (along z axis) and z = 0 (along Vz axis), with widths of 20 km s −1 and 0.3 kpc, respectively.The error bars show the numerical noise.The red lines show the Gaussian profiles fitting to the visually selected intersection regions, with the Gaussian peaks indicated by the light green dashed lines.The intersections, i.e. the Gaussian peaks are marked in the left two columns with the light green dots.

Figure 6 .
Figure 6.Vertical potential derived from the interpolation of the phase snail intersections measured from the test particle simulation.Lower panels: The blue lines show the model potential measured at R = Rg, while the red lines display the potential at R = Rpot that is corrected from the guiding center radius with the relation shown in Fig. 3.The red dots are derived from the interpolation of the intersections between the snails and the z/Vz axes.Upper panels: The relative residuals between the potentials derived from the snail intersections and the model potentials at radii of R = Rg (blue stars) and R = Rpot (red squares).

Figure 7 .
Figure 7. Same as Fig. 5 but for the measurement of the intersections for the observed phase snails.The pixel size is 0.025 kpc × 1.5 km s −1 .

Figure 8 .
Figure 8.The same as Fig. 3 but for the observational data.The bin width of each point is 0.5 kpc.The guiding center radius is calculated according to the rotation curve from Huang et al. (2016).

Figure 9 .
Figure 9.Comparison between the obtained vertical potential of the snail intersections and the popular Milky Way potentials.The red points are the interpolated snail intersections utilizing Eq. 6.The colored lines are the vertical potential profiles from three popular Milky Way potentials, i.e. blue for the Model I in Irrgang et al. (2013), red for the potential from Bovy (2015) (MWPotential2014 in galpy), and green for the potential from McMillan (2017).These vertical potential profiles have been radially shifted according to the different Solar radii (R0) used.The solid lines show the potential at each Rg bin center (as labelled), while the dashed lines display the potential at the corresponding potential radii Rpot (as labelled), corrected from the Rg values according to Fig. 8.apply Gaussian priors for the disc scale lengths, the model predicted values are 4.22 ± 0.27 kpc and 3.99 ± 0.57 kpc, for the thin and thick discs, respectively.As our interpolated snail intersections cover the low vertical regions, our model gives a weak constraint on the thick disc, indicated by a peak of posterior distribution of M d2 close to the lower boundary.Therefore, the overestimation of the vertical potential at outer radii indicates a shallower radial variation, i.e. a larger scale length of the thin disc than its prior.

Figure 10 .
Figure 10.Vertical potential profiles from the global mass modelling with the interpolated snail intersections.The interpolated snail potential data points (red points) are fitted with a model consisting of a double exponential disc and an NFW dark matter halo (black lines), using the corresponding potential radii Rpot of the phase snails.The MWPotential2014 from Bovy (2015) is shown as the red lines for comparison.The solid lines show the potential at each Rg bin center (as labelled), while the dashed lines display the potential at Rpot , corrected from the Rg values using the relation in Fig. 8. (2013); Zhang et al. (2013); Bienaymé et al. (2014); Piffl et al. (2014); Guo et al. (2020); Li & Widrow (2021), as shown in the upper and lower panels of Fig. 14, respectively.The red points and the black solid lines are the same as those in Fig. 12.As shown in the upper panel of Fig.14, our model predicted vertical potential profile, i.e. the median of all potential profiles derived from the posterior parameter sets, shows quite good agreement with the profiles fromLi & Widrow (2021) andHolmberg & Flynn (2000).However, consider-

Figure 11 .
Figure 11.Same as Fig. 7 but for the radial bin at Rg = 7.91 kpc.This phase snail is related to the potential radius of R0.One additional outer wrap is measured compared to Fig. 7.

Figure 12 .
Figure 12.Model predicted local vertical potential (top panel) and total surface density profile (bottom panel) by fitting the interpolated snail intersections (red points).The black lines show the median profiles derived from the posterior parameter sets, while the shaded regions indicate the 1σ errors.The colored lines display the vertical potential at the Solar radius from three popular Milky Way potentials for comparisons.

Figure 13 .
Figure 13.Posterior probability distributions of parameters for the mass modelling of the snail at Rpot = R0.The parameters from the left to right are the stellar surface density (Σ⋆), the stellar disc scale height (Zh) and the local dark matter density (ρdm), respectively.The black dashed lines in the histograms show the medians and 1σ errors of the marginalized 1D distributions, while the red dashed lines indicate the best-fit parameters.
Figure 14.Comparing the model predicted vertical potential (top panel) and the total surface density profile (bottom panel) at Rpot = R0 with previous results in the literature (as labelled).The black lines show our modeling results by fitting the interpolated snail intersections (red points), with shaded regions indicating the 1σ errors.Our results best agree with the profiles from Holmberg & Flynn (2000) (black dashed lines) and Li & Widrow (2021) (blue solid lines).

Figure 15 .
Figure 15.The number density difference between the positive and negative vertical height (left) and vertical velocity (right) at different guiding center radii (as labelled in each row).The black dashed lines are the Gaussian smoothing profiles of the data points.The light green vertical lines indicate the snail intersections measured in Fig. 7.

Table 1 .
The corresponding angular momenta (Lz) and the potential radii (Rpot ) for the Rg bins used in our sample.