# Vertical structure of Galactic disk kinematics from LAMOST K giants

Ping-Jie Ding  CAS Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Beijing 100101, China  Xiang-Xiang Xue CAS Key Laboratory of
Optical Astronomy, National Astronomical Observatories, Beijing 100101, China
Chengqun Yang CAS Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Beijing 100101, China  School of Astronomy and Space Science,
University of Chinese Academy of Sciences, 19A Yuquan Road, Shijingshan
District, Beijing 100049, P.R.China  Shanghai Astronomical Observatory, 80
Nandan Road, Shanghai 200030, People s Republic of China  Gang Zhao CAS Key
Laboratory of Optical Astronomy, National Astronomical Observatories, Beijing
100101, China  School of Astronomy and Space Science, University of Chinese
Academy of Sciences, 19A Yuquan Road, Shijingshan District, Beijing 100049,
P.R.China  Lan Zhang  CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Beijing 100101, China  Zi Zhu  School of Astronomy
and Space Science, Nanjing University, Nanjing 210023, China  Key Laboratory
for Modern Astronomy and Astrophysics (Nanjing University), Ministry of
Education, Nanjing 210023, China

###### Abstract

We examine the vertical structure of Galactic disk kinematics over a
Galactocentric radial distance range of $R=5$-15 kpc and up to 3 kpc away from
the Galactic plane, using the K-type giants surveyed by LAMOST. Based on
robust measurements of three-dimensional velocity moments, a wobbly disk is
detected in a phenomenological sense. An outflow dominates the radial motion
of the inner disk, while in the outer disk there exist alternate outward and
inward flows. The vertical bulk velocities is a combination of breathing and
bending modes. A contraction-like breathing mode with amplitudes increasing
with the distance to the plane and an upward bending mode dominate the
vertical motion outside $R_{0}$, and there are reversed breathing mode and
bending mode at $R<R_{0}$, with amplitudes much smaller than those outside
$R_{0}$. The mean azimuthal velocity decreases with the increasing distance to
the plane, with gradients shallower for larger $R$. Stars in the south disk
are rotating faster than stars in the north. The velocity ellipsoid
orientation differs between different $R$: in the range of
$5<R<9\,\mathrm{kpc}$, the gradient of the tilt angle with respect to
$\arctan(Z/R)$ decreases from $\sim 0.83$ for the inner disk to $\sim 0.52$
for the outer disk; within $9<R<15\,\mathrm{kpc}$, the tilt of velocity
ellipsoid deviates from vertical antisymmetry. A clear flaring signature is
found for both north and south disks based on the observed vertical structures
of velocity ellipsoid.

stars: kinematics and dynamics – Galaxy: kinematics and dynamics – Galaxy:
disk

##  1 Introduction

Our host galaxy, the Milky Way, is a typical disk galaxy. The vertical
structure of the disk kinematics is important for our understanding of the
Galactic formation and evolution. In the classical characterization of the
Galactic disk, the distribution of stellar kinematics along the vertical
distance to the Galactic plane is monotonous and symmetric. Owing to
observations carried out by large surveys in recent years such as the Sloan
Digital Sky Survey (SDSS; York et al. 2000), the RAdial Velocity Experiment
(RAVE; Steinmetz et al. 2006), the Large Aperture Multi-Object Fibre
Spectroscopic Telescope (LAMOST; Cui et al. 2012), and the Gaia mission (Gaia
Collaboration et al. 2016), more and more details in the velocity profiles
have become apparent.

The bulk motions in the Galactocentric radial and vertical directions and
their variations with respect to the vertical height above and below the plane
have been studied using different tracers. Within the disk, there is evidence
of a wobbly radial velocity ($V_{R}$) along the distance to the plane.
Investigating velocity profiles for F-type stars sampled from the LAMOST data
(Zhao et al. 2006; Cui et al. 2012; Zhao et al. 2012) within Galactocentric
radii $7.8<R<9.8\,\mathrm{kpc}$ and $\pm 2\,\mathrm{kpc}$ from the plane,
Carlin et al. (2013) detected a negative mean radial velocity near the plane,
corresponding to an inward radial flow. They also found that the mean of
$V_{R}$ increases with the increasing of $|Z|$. Based on RAVE (Steinmetz et
al. 2006; Siebert et al. 2011) red-clump stars at $|Z|<2\,\mathrm{kpc}$ and
$6<R<10\,\mathrm{kpc}$, Williams et al. (2013) found an outflow with
$V_{R}=8-10\,\mathrm{km\,s}^{-1}$ for $0<Z<1\,\mathrm{kpc}$ and an inward
motion with $V_{R}=-10\,\mathrm{km\,s}^{-1}$ at $R=9\,\mathrm{kpc}$ and
$-1<Z<-0.5\,\mathrm{kpc}$. Wang et al. (2018) presented an analysis of
kinematics of K giant stars selected from the LAMOST catalog within
$|Z|<2\,\mathrm{kpc}$ and found that the $V_{R}$ above the plane is higher
than that below the plane at $R\sim 10-11\,\mathrm{kpc}$. Similar north-south
asymmetries in $V_{R}$ in the outer disk were also detected by Wang et al.
(2019) and Wang et al. (2020) based on the LAMOST red clump stars. Sampling
Gaia-LAMOST AFGK dwarf stars restricted to 1 kpc from the Sun, Ding et al.
(2019) found a $V_{R}$ profile with a gradient $dV_{R}/dZ\sim
29\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ across the plane.

The vertical motion of the disk can be decomposed into a breathing mode motion
and a bending mode motion. The former is defined by a vertical pattern with
odd parity in the $V_{Z}$ distribution with respect to $Z$, and the latter is
a pattern with even parity in the $Z$-$V_{Z}$ plane (Weinberg et al., 1991).
Both the breathing and the bending modes have been detected in the solar
neighborhood, though with differences between works. Widrow et al. (2012)
found that the vertical motion of the SDSS-SEGUE(Aihara et al. 2011; Yanny et
al. 2009) stars resemble that of a breathing mode motion, in which stars above
the plane are on average moving upwards, while those below the plane,
downwards. Carlin et al. (2013) argued that stars both above and below the
plane are moving towards the plane, as a compression movement. A more complex
structure was detected by Williams et al. (2013). They found opposite
breathing modes inside and outside the solar circle: the $V_{Z}$ field implies
a rarefaction at $R<R_{0}$, whereas a contraction at $R>R_{0}$. The inner-disk
breathing mode found by Wang et al. (2018) is the inverse of that found by
Williams et al. (2013). Moreover, they detected a bending motion in the range
of $6<R<13\,\mathrm{kpc}$, which is downward interior to the solar circle and
upward at $R\gtrsim 9\,\mathrm{kpc}$. The findings of Wang et al. (2019) and
Wang et al. (2020) agree with a significant upward bending motion in the outer
disk. Using data from RAVE DR5 (Kunder et al., 2017) and the Tycho-Gaia
astrometric solution catalogue (Michalik et al., 2015), Carrillo et al. (2018)
claimed that a rarefaction-like breathing mode dominates the $V_{Z}$ pattern
inside $R_{0}$, while the pattern of the outer disk is dominated by a downward
bending mode within $|Z|\lesssim 1\,\mathrm{kpc}$. A similar combination of
vertical modes was reported by Carrillo et al. (2019) based on Gaia DR2
astrometric and line-of-sight velocity information, while a positive bending
mode was observed at $R\thickapprox 13\,\mathrm{kpc}$. In the work of López-
Corredoira et al. (2020), the authors found that both the inner and the outer
disks present a combination of breathing and bending mode motions with Gaia
DR2: a strong bending mode exists across the disk, the velocity of which
changes from negative at $R<R_{0}$ to positive at $R>R_{0}$ similar to the
pattern found by Wang et al. (2018); a relatively weak breathing mode shows a
smooth distribution over the disk, with amplitudes highly dependent on the
distance to the plane. Recently, Gaia Collaboration et al. (2021) adopted
proper motions of stars in the direction of the Galactic anticenter from the
Gaia EDR3 (Gaia Collaboration et al. 2020; Lindegren et al. 2020) and detected
upward vertical motions in both north and south disks. They found that the
$Z>0$ stars exhibit faster vertical motions than the $Z<0$ stars by up to
$\sim 2\,\mathrm{km\,s}^{-1}$ at $R<10\,\mathrm{kpc}$, while at
$R>10\,\mathrm{kpc}$ the vertical velocities of $Z<0$ stars are larger than
those of $Z>0$ stars by up to $\sim 6\,\mathrm{km\,s}^{-1}$, implying a
combination of an alternate breathing mode motion and an upward bending mode
motion.

Compared to the behaviors of $V_{R}$ and $V_{Z}$, the vertical structure of
the azimuthal velocities ($V_{\phi}$) is more regular, with a negative
$V_{\phi}$-$|Z|$ trend on both sides of the disk (Smith et al. 2012; Karaali
et al. 2014; Recio-Blanco et al. 2014; Guiglion et al. 2015; Jing et al. 2016;
Carrillo et al. 2018). The vertical gradient in the azimuthal velocities is
dependent on the tracers. For instance, sampling solar neighborhood dwarf
stars in the range of $7<R<9\,\mathrm{kpc}$ and $|Z|<2\,\mathrm{kpc}$, Smith
et al. (2012) found that the gradient of the mean $V_{\phi}$ with respect to
$|Z|$ is around $-15$ to $-40\,\mathrm{km\,s^{-1}\,kpc^{-1}}$. Using 130043
F/G-type dwarf stars within $6.5<R<9.5\,\mathrm{kpc}$ and
$0.1<|Z|<3\,\mathrm{kpc}$ from the LAMOST survey and the SDSS, Jing et al.
(2016) obtained the gradient as $-18.5\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ for the
metallicity range $\rm{[Fe/H]}>-0.1$ and
$-14.2\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ for $-0.8<\rm{[Fe/H]}<-0.6$. In most of
the previous research, the patterns of $V_{\phi}$ were approximately symmetric
about the Galactic plane. However, some small-amplitude asymmetric structures
have been detected in recent studies. Based on proper motions of the red clump
giants from the PPMXL survey (Roeser et al., 2010), López-Corredoira et al.
(2014) found that the rotation speed of the south disk is higher than that in
the north in the range of $8<R<13\,\mathrm{kpc}$, though with large error
bars. Wang et al. (2018) claimed that the K giants within $R\geqslant
10\,\mathrm{kpc}$ and $Z\geqslant 0.5\,\mathrm{kpc}$ exhibit larger rotation
speed in the south disk than in the north, while at $R\thickapprox
8-9\,\mathrm{kpc}$ and $|Z|\leqslant 0.5\,\mathrm{kpc}$ northern stars are
rotating faster than the southern. In the work of Ding et al. (2019), the
authors found that dwarfs at $-0.5<Z<0\,\mathrm{kpc}$ are rotating slightly
faster than those at $0<Z<0.5\,\mathrm{kpc}$. The faster-rotating of the
southern stars was also evidenced by Gaia Collaboration et al. (2021). They
claimed that the rotation of $Z<0$ stars leads that of $Z\geq 0$ stars at
$R\gtrsim 11$ kpc, typically by up to $10\,\mathrm{km\,s}^{-1}$.

In addition to the bulk motions, the velocity ellipsoid as a function of the
vertical height provides information on the structure and evolution process of
the disk. The velocity dispersion is known to be caused by the disk heating
due to fluctuations in the gravitational field. An increment with the distance
to the plane of all the velocity dispersions (i.e. $\sigma_{R}$, $\sigma_{Z}$
and $\sigma_{\phi}$) has been recognized in earlier works (e.g. Soubiran et
al. 2003; Moni Bidin et al. 2012; Smith et al. 2012; Karaali et al. 2014; Ding
et al. 2019; Mackereth et al. 2019), indicating that the heating process is
stronger for higher disk stars. Moreover, it has been found that the
covariance between $V_{R}$ and $V_{Z}$, often parameterized by the tilt angle
($\alpha_{\rm{tilt}}$) of the orientation of the velocity ellipsoid, is
antisymmetric about the plane, and the orientation of the velocity ellipsoid
stays between an horizontal one ($\alpha_{\rm{tilt}}\equiv 0$) and a spherical
one ($\alpha_{\rm{tilt}}=\arctan(Z/R)$). In the work of Binney (2014), the
authors found $\alpha_{\rm{tilt}}\sim 0.8\arctan(Z/R)$ by using RAVE stars
lying within $\sim 2\,\mathrm{kpc}$ of the Sun, suggesting a near-spherical
orientation of the velocity ellipsoid. The findings of Büdenbender et al.
(2015) and Ding et al. (2019) based on nearby dwarf stars, which gave
$\alpha_{\rm{tilt}}=(0.90\pm 0.04)\arctan(Z/R_{0})-(0.01\pm 0.005)$ and
$\alpha_{\rm{tilt}}=(1.11\pm 0.11)\arctan(Z/R_{0})-(0.0069\pm 0.0034)$
respectively, agree with a spherical alignment of the tilt in the solar
vicinity. Using data from the Gaia RVS catalog (Cropper et al., 2019) within 5
kpc of the Sun, Everall et al. (2019) found the gradient of
$\alpha_{\rm{tilt}}$ with $\arctan(Z/R)$ is around 0.95. They also claimed
that the tilt is closer to a spherical orientation at larger $R$ and higher
$|Z|$. To the almost contrary, Hagen et al. (2019) found a (near) spherical
orientation of the velocity ellipsoid at $R\lesssim 7\,\mathrm{kpc}$ and a
more cylindrical tilt beyond $R\gtrsim 9\,\mathrm{kpc}$ by sampling stars
within $4\lesssim R\lesssim 13\,\mathrm{kpc}$ and $|Z|\lesssim
3.5\,\mathrm{kpc}$ from Gaia DR2.

The structural complexity in the vertical distribution of stellar velocities
and the inconsistency in the earlier findings call for a reanalysis of
vertical profiles of velocity moments in a phenomenological sense, by using
favorable tracers obtained from developed observations. In this work, we focus
on how the measured three-dimensional velocity moments modify the vertical
patterns of disk kinematics, by sampling K giant stars with fundamental
stellar parameters and line-of-sight velocities from the recent released
LAMOST catalog. With the help of the high-precision measurements of proper
motions provided by the Gaia EDR3, we examine in detail the variation of
velocity moments with the distance to the plane in a large volume around the
Sun and analyze the Galactic structure inferred from the observed kinematics.

In Section 2, we introduce our sample of K giant stars, for which we have
LAMOST line-of-sight velocities, distances, and metallicities, matched to Gaia
EDR3 proper motions. Section 3 gives the vertical distribution of velocity
moments. The difference between the inner and outer disks as well as the
dependence of stellar kinematics on the metallicity will also be presented. In
Section 4, we provide a measurement of the flaring strength of the stellar
disk by using the results in Section 3. We draw our discussion in Section 5
and present a summary of our results in Section 6.

##  2 Sample selection

In this work, we focus on the characterisation of the vertical structure of
velocity moments across the Galactic disk. The K-type giants are a kind of
long-lived objects with high luminosities. The large sample size and the wide
spatial coverage of the observed K giants make them feasible tracers for the
Galactic kinematics. We sample the K giants from the catalog of LAMOST DR8,
which contains observations from October 24, 2011 to May 27, 2020. The K
giants are selected based on the surface gravities ($\rm{log}\emph{g}$) and
effective temperatures ($T_{\rm{eff}}$) measured by the LAMOST team. The
selection criteria follow the conditions given by Liu et al. (2014) such that
$4000\,\rm{K}<T_{\rm{eff}}<4600\,\rm{K}$ with $\rm{log}\emph{g}<3.5$ and
$4600\,\rm{K}<T_{\rm{eff}}<5600\,\rm{K}$ with $\rm{log}\emph{g}<4$. We obtain
the metallicities ($\rm{[Fe/H]}$) and the line-of-sight velocities
($v_{\rm{los}}$) from the LAMOST data. We reject stars with
$\rm{[Fe/H]}\leqslant-1$ to exclude likely halo stars. It has been noted that
there exists a systematic offset in the measured $v_{\rm{los}}$ of around
$5-7\mathrm{km\,s}^{-1}$ (Gao et al. 2015; Luo et al. 2015; Tian et al. 2015;
Jing et al. 2016; Schönrich & Aumer 2017), which needs to be added to
$v_{\rm{los}}$. We select stars from the LAMOST sample that also have
measurements of line-of-sight velocities in the Gaia data to determine the
offset. The Gaia EDR3 contains line-of-sight velocities for more than 7
million sources (Cropper et al. 2019, Gaia Collaboration et al. 2020), with
systematic bias at a level of $0.1\,\mathrm{km\,s}^{-1}$ (Katz et al., 2019).
We use the Gaia line-of-sight velocities satisfying
$\rm{RV_{-}NB_{-}TRANSITS}>5$ as calibrators, and find that the offset in the
LAMOST measurement is $5.34\,\mathrm{km\,s}^{-1}$ (see Fig. 1).

Figure 1: Line-of-sight velocities of the LAMOST DR8 K giants versus those
provided by the Gaia EDR3. The red line denotes the linear fit with a slope of
1.0.

Figure 2: Distances derived by the method of Carlin et al. (2015) versus the
inverse of the Gaia parallaxes for the $\sim 355000$ K giants with high-
precision trigonometric parallaxes (i.e. $\sigma_{\varpi}/\varpi<0.1$). The
red line denotes the linear fit with a slope of 0.842.

Figure 3: Error histograms for space velocities $V_{R}$, $V_{\phi}$ and
$V_{Z}$ respectively for all the selected K giants.

Figure 4: Velocity uncertainties as functions of $Z$ (left panels) and $R$
(right panels) respectively. The lines mark the medians of the uncertainties
of $V_{R}$ (top panels), $V_{\phi}$ (middle panels) and $V_{Z}$ (bottom
panels) respectively. The shaded regions, limited by the 16 and 84
percentiles, show areas enclosing 68 percent of the stars.

Figure 5: Histograms for the metallicity $\rm{[Fe/H]}$ in the north ($Z>0$)
and south ($Z<0$) disks respectively.

In order to obtain a robust characterisation of stellar motions, we need
information of proper motions and distances. The proper motions are obtained
from the Gaia EDR3 (Gaia Collaboration et al. 2020; Lindegren et al. 2020),
which provides the best astrometric data nowadays. We adopt the criteria
$\rm{ASTROMETRIC_{-}GOF_{-}AL}<3.0$ and
$\rm{ASTROMETRIC_{-}EXCESS_{-}NOISE_{-}SIG}\leqslant 2.0$ to select
astrometrically well-behaved stars with good-fit statistic of the astrometric
solution, coupled with a selection of the Renormalised Unit Weight Error
(RUWE) $<1.4$ (Lindegren et al., 2020) to favour stars with good astrometric
quality. The distances, for the geometry, can be simply estimated by the
inverse of the trigonometric parallax. However, such a method is feasible only
for parallaxes with very high precision (Bailer-Jones et al. 2015;
Astraatmadja & Bailer-Jones 2016). Therefore, we employ the photometric
distances derived from the measurements of distance modulus. To do this, we
adopt the absolute magnitudes estimated using the method of Carlin et al.
(2015), which is tailored specifically to the spectroscopically derived
stellar parameters of LAMOST. Comparison with the apparent magnitude, which is
obtained from the 2MASS point source catalog (Cutri et al. 2003; Skrutskie et
al. 2006) in the $K_{s}$ band (2.16 $\mu$m), yields the distance modulus and
then the distance. In order to correct for the extinction, we adopt
$R(K_{s})=0.306$ from Yuan et al. (2013) as the reddening relative to $E(B-V)$
obtained from Schlegel et al. (1998). We obtain 545227 K giants with
measurements of proper motions and distances.

We compare the derived photometric distances to the inverse of the Gaia
parallax. For the comparison, we select $\sim 355000$ Gaia stars with small
relative parallax uncertainties $\sigma_{\varpi}/\varpi<0.1$ since using
$1/\varpi$ as the distance estimate only works for stars with high precision.
The comparison is shown in Fig. 2, in which the zeropoint of the Gaia
parallax, i.e., 0.017 mas (Lindegren et al., 2020), has been added to the
parallaxes. There is close agreement between the two distance scales, with a
ratio of approximately 0.842 of the photometric distances to the inverse of
parallaxes. In this case, we calibrate the photometric distances to more
proper ones by dividing them by 0.842.

Throughout this paper we use Galactocentric cylindrical coordinates
$(R,\,\phi,\,Z)$ with $R$ pointing outwards from the Galactic center (GC), the
azimuthal $\phi$ in the direction of Galactic rotation, and $Z$ positive
towards the north Galactic pole. The Galactocentric distance of the Sun is
fixed at $R_{0}=8.122\,\mathrm{kpc}$ Gravity Collaboration et al. (2018), and
the vertical offset of the Sun from the Galactic mid-plane is assumed to be
$Z_{\odot}=20.8$ pc (Bennett & Bovy, 2019). The sample is located within
$5<R<15\,\mathrm{kpc}$ and $|Z|<3\,\mathrm{kpc}$. Since the reference of the
Galactocentric system is fixed at the GC, we need to know the velocity of the
Sun. The azimuthal component of the solar velocity ($v_{\odot}$) can be
obtained directly from the proper motion of Sgr A* along the Galactic
longitude, i.e., $\mu_{\rm{SgrA*}}=-6.411\,\mathrm{mas\,yr}^{-1}$ (Reid &
Brunthaler, 2020), which yields $v_{\odot}=246.8\,\mathrm{km\,s}^{-1}$. The
radial and vertical velocities of the Sun are taken from the widely used
$U_{\odot}=11.1\,\mathrm{km\,s}^{-1}$ and
$W_{\odot}=7.25\,\mathrm{km\,s}^{-1}$ given by Schönrich et al. (2010). In our
sample, the velocity is confined to $|V_{R}|<200\,\mathrm{km\,s}^{-1}$,
$-80<V_{\phi}<320\,\mathrm{km\,s}^{-1}$ and $|V_{Z}|<150\,\mathrm{km\,s}^{-1}$
by using the $3\sigma$ criterion to avoid outlier velocities. Finally, a total
number of $N=488209$ K giants remain.

The error histograms for the three-dimensional space velocities of our sample
shown in Fig. 3 indicate good measurement qualities. The medians of $\delta
V_{R}$, $\delta V_{\phi}$ and $\delta V_{Z}$ for all the selected K giants are
5.24, 5.65 and $3.97\,\mathrm{km\,s}^{-1}$ respectively. Figure 4 gives the
velocity uncertainties as functions of vertical distance and Galactocentic
radius respectively. The measurement quality towards the Galactic anticenter
is better than that towards the GC, which means that we can detailed detect
stellar kinematics in the range of $|Z|<3\,\mathrm{kpc}$ from the inner disk
near the solar circle to the distant outer disk.

##  3 Results

First and foremost, we construct the vertical pattern of stellar kinematics
based on all the selected K giants as a function of Galactocentric radius.
Then we compare the kinematics traced by populations with different
metallicities. Figure 5 shows the histograms for the $\rm{[Fe/H]}$ above and
below the plane separately. When investigating the metallicity dependence of
stellar motions, we divide the sample into 5 metallicity bins: the most metal-
poor and metal-rich bins are $-1<\rm{[Fe/H]}<-0.6$ and $0<\rm{[Fe/H]}<0.5$
respectively; stars with intermediate metallicities are equally divided into
bins with a width of 0.2 dex.

###  3.1 Velocity moments in the radial and vertical directions

The radial and vertical components of velocity moments, which carry certain
kinematic details along the distance to the plane, play important roles in
constraining the vertical structure of the disk. We use vertical bins of 0.2
kpc with the step of 0.1 kpc to investigate the velocity moments as functions
of $Z$. Each bin contains no less than 30 stars. When making a comparison
between different $\rm{[Fe/H]}$, the bin width in $Z$ is 1 kpc and the step is
0.5 kpc. In each bin, we obtain the mean and dispersion of velocities in
radial and vertical directions based on the Gaussian distribution (see
Appendix).

####  3.1.1 Radial bulk motions

Figure 6 shows the mean of $V_{R}$, $\langle V_{R}\rangle$, as functions of
$Z$ for populations at different $R$. For stars near the plane, the $\langle
V_{R}\rangle$ changes from positive to negative when we move across the solar
circle and back to positive beyond $R=10\,\mathrm{kpc}$. At higher disk there
are north-south asymmetries in the $\langle V_{R}\rangle$ distribution,
especially in the range of $9<R<13\,\mathrm{kpc}$ where the $\langle
V_{R}\rangle$-$Z$ pattern relates to an outward flow with velocities up to
around $20\,\mathrm{km\,s}^{-1}$ at $Z<0$ and alternate outward and inward
flows with velocities $\lesssim 10\,\mathrm{km\,s}^{-1}$ at $Z>0$. The
$\langle V_{R}\rangle$ increases with the increasing of $|Z|$ within $-2<Z<0$,
with a gradient of $\sim 4-7\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ at
$5<R<11\,\mathrm{kpc}$ and $\sim 0-5\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ at
$11<R<15\,\mathrm{kpc}$. In the north disk, the $\langle V_{R}\rangle$-$Z$
gradient for $R<10\,\mathrm{kpc}$ is positive, shallowing as $R$ increases,
while for $10<R<13\,\mathrm{kpc}$ the gradient stays positive only between
$z\sim 0$ and $0.5\,\mathrm{kpc}$, and reverses beyond $Z=0.5\,\mathrm{kpc}$.
The $Z$-dependence of $\langle V_{R}\rangle$ is nearly vanished for the radial
slice farthest from the GC.

Figure 6: Mean radial velocities as functions of the vertical distance to the
Galactic mid-plane, coupled with the measurement errors. The horizontal axis
is the mean vertical distance of each vertical bin. The range of the
Galactocentric radius, in unit of kpc, is labeled at the top of each panel.

Figure 7: Mean radial velocities as functions of the vertical distance to the
Galactic mid-plane for different metallicities, coupled with the measurement
errors. The black, red, orange, olive and blue dots denote the metallicity
range $-1<\rm{[Fe/H]}<-0.6$, $-0.6<\rm{[Fe/H]}<-0.4$, $-0.4<\rm{[Fe/H]}<-0.2$,
$-0.2<\rm{[Fe/H]}<0$ and $0<\rm{[Fe/H]}<0.5$, respectively. The horizontal
axis is the mean vertical distance of each vertical bin. The range of the
Galactocentric radius, in unit of kpc, is labeled at the top of each panel.

Figure 7 presents the comparison between different $\rm{[Fe/H]}$. There is no
significant relation between $\langle V_{R}\rangle$ and $\rm{[Fe/H]}$ in the
inner disk. In the outer disk within $8<R<12\,\mathrm{kpc}$, metal-poor
populations have higher $\langle V_{R}\rangle$ than metal-rich ones, though
the vertical behaviors of $\langle V_{R}\rangle$ for different $\rm{[Fe/H]}$
are qualitatively consistent. The gap of $\langle V_{R}\rangle$ between
metallicities tends to narrow when $|Z|\gtrsim 1\,\mathrm{kpc}$ and
$R>12\,\mathrm{kpc}$, indicating a mixture of populations at high disk and
large Galactocentric radius.

Figure 8: As in Fig. 6, but for the mean vertical velocities.

Figure 9: Breathing velocities as functions of the height from the Galactic
mid-plane. The range of the Galactocentric radius, in unit of kpc, is labeled
at the top of each panel.

Figure 10: As in Fig. 9, but for the bending velocities.

Figure 11: Breathing velocities as functions of the height from the Galactic
mid-plane for different metallicities. The black, red, orange, olive and blue
dots denote the metallicity range $-1<\rm{[Fe/H]}<-0.6$,
$-0.6<\rm{[Fe/H]}<-0.4$, $-0.4<\rm{[Fe/H]}<-0.2$, $-0.2<\rm{[Fe/H]}<0$ and
$0<\rm{[Fe/H]}<0.5$, respectively. The range of the Galactocentric radius, in
unit of kpc, is labeled at the top of each panel.

Figure 12: As in Fig. 11, but for the bending velocities.

Figure 13: As Fig. 6, but for the radial (black squares) and vertical (red
triangles) velocity dispersions, associated with the ratio of vertical to
radial dispersions (blue circles).

Figure 14: As Fig. 7, but for the radial velocity dispersions.

Figure 15: As Fig. 7, but for the vertical velocity dispersions.

####  3.1.2 Breathing and bending motions

The vertical patterns of the mean vertical velocities $\langle V_{Z}\rangle$
given by Fig. 8 show a north-south asymmetry, which is most notable in the
range $8<R<12\,\mathrm{kpc}$. Figures 9 and 10 give the breathing and bending
velocities as functions of $|Z|$ respectively, which are calculated based on
the the definition in Gaia Collaboration et al. (2018), i.e.,

$\displaystyle V_{\rm{breath}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\langle V_{Z}\rangle^{\rm{North}}-\langle
V_{Z}\rangle^{\rm{South}}),$ $\displaystyle V_{\rm{bend}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\langle V_{Z}\rangle^{\rm{North}}+\langle
V_{Z}\rangle^{\rm{South}}).$ (1)

In the inner disk, there is no significant breathing motion at $|Z|\lesssim
2\,\mathrm{kpc}$; beyond 2 kpc we find positive $V_{\rm{breath}}$ as large as
around $5\,\mathrm{km\,s}^{-1}$, suggesting a rarefaction effect at high disk.
The phase of the breathing mode reverses outside the solar circle: the
vertical motion is dominated by a contraction movement, whose amplitude shows
dependence on both $|Z|$ and $R$: at $|Z|\sim 1\,\mathrm{kpc}$, the
$|V_{\rm{breath}}|$ is $\sim 5\,\mathrm{km\,s}^{-1}$ within
$8<R<10\,\mathrm{kpc}$ while stays $\sim 0-2\,\mathrm{km\,s}^{-1}$ within
$10<R<15\,\mathrm{kpc}$; at $|Z|\sim 2\,\mathrm{kpc}$, the $|V_{\rm{breath}}|$
increases with $R$ when $R<11\,\mathrm{kpc}$ while the trend reverses beyond
$R=11\,\mathrm{kpc}$. In the range between $R=9$ and 15 kpc, the amplitude of
the breathing motion increases towards higher disk, the gradient of
$|V_{\rm{breath}}|$ with $|Z|$ decreasing from $\sim
3.7\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ at $9<R<10\,\mathrm{kpc}$ to $\sim
1.4\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ at $13<R<15\,\mathrm{kpc}$.

Similar to the breathing velocity, the sign of the bending velocity reverses
when we move across the solar circle (see Fig. 10). In the inner disk, the
bending motion is trivial. The mean of $V_{\rm{bend}}$ is around
$-1\,\mathrm{km\,s}^{-1}$ at $R<8\,\mathrm{kpc}$. In the outer disk, the
$V_{\rm{bend}}$ stays positive (up to $\sim 5\,\mathrm{km\,s}^{-1}$) within
$|Z|<2.5\,\mathrm{kpc}$, suggesting a bending movement towards the north.
There is a flat or positive $V_{\rm{bend}}$-$|Z|$ gradient at
$|Z|<1\,\mathrm{kpc}$. When we move to $|Z|>1\,\mathrm{kpc}$, a negative
$V_{\rm{bend}}$-$|Z|$ gradient dominates the bending mode pattern.

Figures 11-12 present the breathing and bending velocities for populations
with different metallicities. The differences between the amplitudes of the
breathing motions for different $\rm{[Fe/H]}$ are only evident in the solar
vicinity, where the $|V_{\rm{breath}}|$ increases with the increasing of
$\rm{[Fe/H]}$. In the region far away from the Sun, the $|V_{\rm{breath}}|$
for different $\rm{[Fe/H]}$ are nearly consistent. As for the bending motion,
we find a positive $V_{\rm{bend}}$-$\rm{[Fe/H]}$ trend at
$8<R<9\,\mathrm{kpc}$, similar to the metallicity trend of
$|V_{\rm{breath}}|$. Outside of this slice, the $V_{\rm{bend}}$ is almost
independent on $\rm{[Fe/H]}$, except for that in the range of
$7<R<12\,\mathrm{kpc}$ the $V_{\rm{bend}}$ of the most metal-poor population
is lower than that of other populations.

####  3.1.3 Velocity dispersions in the radial and vertical directions

The radial and vertical velocity dispersions are presented in Fig. 13. Both
$\sigma_{R}$ and $\sigma_{Z}$ increase with the increasing distance to the
plane, which means that populations further away from the plane have been
stronger affected by the heating mechanisms so that been diffused through a
larger fraction of phase space. The radial dispersions on both sides of the
disk are approximately comparable, only with small differences at
$R<10\,\mathrm{kpc}$ and $|Z|<2\,\mathrm{kpc}$ that the $\sigma_{R}$ above the
plane is on average larger than that below the plane by no more than a few
$\mathrm{km\,s}^{-1}$. A more significant north-south asymmetry exists in
$\sigma_{Z}$ within $7<R<13\,\mathrm{kpc}$ such that $\sigma_{Z}$ of $Z<0$
stars are larger than that of $Z>0$ stars by up to $\sim
15\,\mathrm{km\,s}^{-1}$.

The velocity dispersions decrease towards the outer disk. In comparison, there
is no pronounced trend between $\sigma_{Z}/\sigma_{R}$ and $R$ (see Fig. 13).
The mean $\sigma_{Z}/\sigma_{R}$ is 0.49 at $Z\sim 0$, then increases with the
increasing $|Z|$, notably when $|Z|\lesssim 1.5\,\mathrm{kpc}$. The mean
$\sigma_{Z}/\sigma_{R}$ is 0.73 at $Z<-1.5\,\mathrm{kpc}$ and 0.67 at
$Z>1.5\,\mathrm{kpc}$.

The metallicity dependence of $\sigma_{R}$ and $\sigma_{Z}$ is shown in Figs.
14 and 15 respectively. The gradient of the dispersions with respect to $|Z|$
is larger for populations with higher metallicities. The increasing of both
$\sigma_{R}$ and $\sigma_{Z}$ with the decreasing $\rm{[Fe/H]}$ is notable
within $5<R<10\,\mathrm{kpc}$, while when we move towards the Galactic
anticenter, the dependence of dispersions on $\rm{[Fe/H]}$ weakens, especially
for the radial component.

####  3.1.4 Tilt of the velocity ellipsoid

We present the tilt angles of the velocity ellipsoid in Fig. 16, which are
calculated by the covariance between $V_{R}$ and $V_{Z}$ via
$\alpha_{\rm{tilt}}\equiv\frac{1}{2}\arctan[2\sigma_{RZ}/(\sigma_{R}^{2}-\sigma_{Z}^{2})]$.
The panels of $5<R<9\,\mathrm{kpc}$ show a clear positive trend between
$\alpha_{\rm{tilt}}$ and $Z$. The vertical patterns of $\alpha_{\rm{tilt}}$
can be fitted well by $\alpha_{\rm{tilt}}=(0.827\pm
0.040)\arctan(Z/R)-(0.007\pm 0.011)$, $\alpha_{\rm{tilt}}=(0.793\pm
0.034)\arctan(Z/R)-(0.0032\pm 0.0061)$ and $\alpha_{\rm{tilt}}=(0.520\pm
0.041)\arctan(Z/R)+(0.0370\pm 0.0064)$ at $5<R<7\,\mathrm{kpc}$,
$7<R<8\,\mathrm{kpc}$ and $8<R<9\,\mathrm{kpc}$, respectively, indicating that
the orientation of the velocity ellipsoid is near-spherical in the inner disk
and turns more horizontal when we move across the solar circle.

Figure 16: As Fig. 6, but for the tilt angles.

In the range of $9<R<15\,\mathrm{kpc}$, the
$\alpha_{\rm{tilt}}$-$\arctan(Z/R)$ gradient ranges between $\sim 0.4$ and
$\sim 1.2$ for $-1<Z<3\,\mathrm{kpc}$. It is interesting to find that the
vertical behavior of $\alpha_{\rm{tilt}}$ reverses at
$Z\lesssim-1\,\mathrm{kpc}$, which breaks the antisymmetry in the
$\alpha_{\rm{tilt}}$-$Z$ plane. Figure 17 shows the tilt angles for different
metallicities. The vertical patterns of $\alpha_{\rm{tilt}}$ for different
$\rm{[Fe/H]}$ are basically consistent.

Figure 17: As Fig. 7, but for the tilt angles.

###  3.2 Vertical structures of azimuthal velocities

In this subsection, we focus on the azimuthal velocities as functions of
distance to the plane. We take a bin of 0.5 kpc and a step of 0.25 kpc to
construct the vertical structure of velocity moments in the azimuthal
direction. Each bin contains no less than 100 stars. We employ a non-Gaussian
distribution function to fit the observed distribution of the azimuthal
velocity (see Appendix) and obtain the first and second moments from the
distribution function.

Figure 18: As Fig. 6, but for the mean azimuthal velocities.

Figure 19: As Fig. 7, but for the mean azimuthal velocities.

Figure 18 shows the mean azimuthal velocity, $\langle V_{\phi}\rangle$, as
functions of $Z$ in different radial slices. The $\langle V_{\phi}\rangle$ at
$Z\sim 0$ shows no significant radial trend, suggesting a near-flat rotation
curve in the plane. The $\langle V_{\phi}\rangle$ decreases as we move away
from the plane. Table 1 gives the absolute values of the vertical gradient in
$\langle V_{\phi}\rangle$ above and and below the plane respectively. The
$|d\langle V_{\phi}\rangle/dZ|$ in the south disk is smaller than that in the
north, especially in the range of $5<R<12\,\mathrm{kpc}$, suggesting that
stars below the plane are on average rotating faster than stars above. Figure
19 presents the $\langle V_{\phi}\rangle$ patterns for different populations.
Generally, the azimuthal velocities of populations with $-1<\rm{[Fe/H]}<-0.4$
increase towards the outer disk, while the trend reverses for
$-0.2<\rm{[Fe/H]}<0.5$.

Figure 20: As Fig. 6, but for the azimuthal velocity dispersions.

Figure 21: As Fig. 7, but for the azimuthal velocity dispersions. Table 1: Absolute values of the vertical gradient in the mean azimuthal velocities in the north ($Z>0$) and south ($Z<0$) disks respectively within different radius slices. $R\,(\mathrm{kpc})$ |  | $(5,\,7)$ | $(7,\,8)$ | $(8,\,9)$ | $(9,\,10)$ | $(10,\,11)$ | $(11,\,12)$ | $(12,\,13)$ | $(13,\,15)$
---|---|---|---|---|---|---|---|---|---
$|d\langle V_{\phi}\rangle/dZ|\,(\mathrm{km\,s^{-1}\,kpc^{-1}})$ | $Z>0$ | $28.21\pm 0.79$ | $27.64\pm 0.57$ | $24.9\pm 1.0$ | $21.59\pm 0.79$ | $20.35\pm 0.73$ | $16.45\pm 0.39$ | $11.00\pm 0.29$ | $4.92\pm 0.30$
| $Z<0$ | $26.79\pm 0.84$ | $24.74\pm 0.74$ | $18.6\pm 1.2$ | $14.4\pm 1.3$ | $14.49\pm 0.40$ | $14.73\pm 0.98$ | $10.4\pm 1.0$ | $5.21\pm 0.98$

The vertical patterns of the azimuthal dispersion $\sigma_{\phi}$ and their
dependence on $\rm{[Fe/H]}$ are given by Figs. 19 and 20 respectively. As
expected, the $\sigma_{\phi}$ increases towards higher disk and lower
metallicity as a whole. The nouth-south asymmetry in $\sigma_{\phi}$ is less
pronounced than that in $\langle V_{\phi}\rangle$. In the innermost slice, the
$\sigma_{\phi}$ is asymmetric within $|Z|\lesssim 1\,\mathrm{kpc}$, with
larger values for north disk stars. The reverse applies for the
$\sigma_{\phi}$ pattern of the the outer disk.

##  4 An implication of the vertical structure of stellar kinematics: the
flare of the Galactic disk

The vertical structure of stellar kinematics gives us an insight into the
stellar distribution in the vertical direction. In this section, we use the
observed velocity ellipsoid shown in Sect. 3 to find the disk scaleheights
$h_{z}$ at different $R$, which allows us to detect the flaring strength of
the disk.

We employ the method from Moni Bidin et al. (2012) and López-Corredoira et al.
(2020) to estimate the $h_{z}$ and its radial behavior. The basic idea is to
fit the expression for the surface density $\Sigma(Z)$ derived from the Jeans
equations to the measured velocity dispersions and covariance:

$\displaystyle\Sigma(Z)$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
G}\Big{[}k_{1}\cdot\int_{0}^{Z}\sigma_{R}^{2}\rm{d}Z+k_{2}\cdot\int_{0}^{Z}\sigma_{\phi}^{2}\rm{d}Z+k_{3}\cdot\sigma_{RZ}$
(2)
$\displaystyle+\frac{\sigma_{Z}^{2}}{h_{z}}-\frac{\partial\sigma_{Z}^{2}}{\partial
Z}+|Z|\sigma_{RZ}\frac{\partial}{\partial
R}\Big{(}\frac{1}{h_{z}}\Big{)}+\int_{0}^{Z}|Z|\frac{\sigma_{R}^{2}}{R}\frac{\partial}{\partial
R}\Big{(}\frac{1}{h_{z}}\Big{)}\rm{d}Z$
$\displaystyle+\int_{0}^{Z}|Z|\sigma_{R}^{2}\frac{\partial^{2}}{\partial
R^{2}}\Big{(}\frac{1}{h_{z}}\Big{)}\rm{d}Z\Big{]}$ $\displaystyle=$
$\displaystyle 2\rho(R,\,Z=0)\cdot
h_{z}(R)\cdot\big{[}1-e^{-Z/h_{z}(R)}\big{]}.$

The constants $k_{1}$, $k_{2}$ and $k_{3}$ are defined by

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle\frac{3}{R_{0}\cdot
h_{R}}-\frac{2}{h_{R}^{2}},$ $\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{R_{0}\cdot h_{R}},$ $\displaystyle k_{3}$
$\displaystyle=$ $\displaystyle\frac{3}{h_{R}}-\frac{2}{R_{0}},$ (3)

where $h_{R}$ is the radial scalelength. In this work we fix $h_{R}\equiv
4.9\,\mathrm{kpc}$ from the result of López-Corredoira et al. (2020).

We derive the fitting values for $h_{z}$ using an iteration method. Using the
initial values of $h_{z}(R)$, the corresponding derivatives of $1/h_{z}$,
namely $\frac{\partial}{\partial R}\Big{(}\frac{1}{h_{z}}\Big{)}$ and
$\frac{\partial^{2}}{\partial R^{2}}\Big{(}\frac{1}{h_{z}}\Big{)}$, and the
$\Sigma(Z)$ can be derived for every radial slice. Then the new values of
$h_{z}(R)$ are determined from Eq. (2) using the least square method. We adopt
the expression $h_{z}(R)=0.533+0.103(R-R_{0})$ from López-Corredoira et al.
(2020) to calculate the initial values for $h_{z}$. The iteration procedure is
repeated until convergence is achieved.

When using Eq. (2) to estimate $h_{z}$ we adopt two assumptions behind the
Jeans equations, namely the stationary of the disk and the axisymmetry of the
system. There is evidence that some non-axisymmetries exist in the disk, which
means that the Jeans equations do not strictly describe the Galactic dynamics
and the $h_{z}(R)$ derived from Eq. (2) is an approximation of the scaleheight
for the observation range.

Figure 22 shows the estimated $h_{z}$ derived from the measured velocity
ellipsoids in the north and south disks respectively. (For $Z<0$ we use $|Z|$
instead of $Z$ and $-\sigma_{RZ}$ instead of $\sigma_{RZ}$ in Eq. (2).) It is
interesting to find that the $h_{z}$ traced by south disk stars is larger than
that traced by north disk stars, with difference decreasing with the
increasing $R$. Nonetheless, the flaring features in the north and south disks
are similar, in which the noticeable flaring begins at $R\sim
8-9\,\mathrm{kpc}$. We use a second-order polynomial fit to constrain the
flaring. The new expressions for $h_{z}(R)$ are $h_{z}=(0.265\pm
0.031)+(0.100\pm 0.025)(R-R_{0})+(0.0261\pm
0.0076)(R-R_{0})^{2}\,\mathrm{kpc}$ for the north disk, and $h_{z}=(0.471\pm
0.029)+(0.091\pm 0.023)(R-R_{0})+(0.0186\pm
0.0069)(R-R_{0})^{2}\,\mathrm{kpc}$ for the south disk.

Figure 22: Scale heights as functions of the Galactocentric radius in the
north (black squares) and south (blue triangles) disks respectively. The fits
of the scale heights are marked by the black and blue lines respectively.

##  5 Discussion

Due to high luminosities and a large sample size, the K giant stars are
increasingly being used as tracers to probe Galactic kinematics and dynamics
over a large section of the Milky Way. The huge number spectra observed by the
LAMOST survey yield an important K giants sample with stellar parameters and
distance estimates. With the help of the high-precision astrometric data
provided by the Gaia EDR3, we can characterize the stellar disk in a wider
range of vertical distance with unprecedented precision. We find north-south
asymmetries in all the three components of mean velocities, where the vertical
patterns of the radial and vertical components are more complex than the
azimuthal components. The non-zero $\langle V_{R}\rangle$ and $\langle
V_{Z}\rangle$ and their variations with $Z$ in different radial slices suggest
that there exist perturbations with internal or external origins acting on
stellar motions.

We find non-zero mean radial velocities ranging from
$\sim-10\,\mathrm{km\,s}^{-1}$ to $\sim 20\,\mathrm{km\,s}^{-1}$, which are
dependent on the distance to the plane. The $\langle V_{R}\rangle$ decreases
towards larger $R$ as a whole, coupled with shallowing vertical gradients.
Considering different metallicity bins, there is evidence that the $\langle
V_{R}\rangle$ of metal-rich stars decreases with $R$ more steeply than that of
metal-poor stars in the range of $8<R<12\,\mathrm{kpc}$, which is in agreement
with the finding of Wojno et al. (2018) from the RAVE FGK turn-off stars. The
observed radial flow and its alternate direction near the plane are likely
generated by the wave-like spiral structures in the plane (Faure et al. 2014;
Monari et al. 2016). Faure et al. (2014) argued that the effect of the spirals
on the stellar radial velocities keeps strong up to $Z\sim\pm
0.5\,\mathrm{kpc}$. Besides, Monari et al. (2014) found that the central bar
could produce radial bulk motions within $R<10\,\mathrm{kpc}$ and give rise to
a negative radial gradient in the radial bulk velocities in the range from
$Z=0$ to $Z\sim 2\,\mathrm{kpc}$. Within $8<R<12\,\mathrm{kpc}$, populations
with different $\rm{[Fe/H]}$ respond differently to the perturbations, with
the tendency for the metal-rich populations to show a larger fractional
contribution to the disturber than the metal-poor ones, though all of them
have consistent vertical trends. The behaviors of $\langle V_{R}\rangle$ agree
with the hypothesis that the in-plane substructures, mainly including the
spirals and the bar, play roles in producing the structure of radial bulk
velocities in the disk. Nonetheless, we cannot exclude the possibility that an
external perturber, such as the passage of a satellite, causes the wobbly
radial bulk motions.

Unlike the radial component of stellar bulk motion, the vertical component,
which gets translated as a combination of breathing and bending mode motions,
is stronger outside $R_{0}$ than inside. The breathing motion detected in the
outer disk confirms the findings of Williams et al. (2013) that there is a
notable vertical contraction at $R>R_{0}$. The inner-disk rarefaction-like
vertical motion, which was reported by some of the earlier works (e.g.
Williams et al. 2013; Carrillo et al. 2018; Carrillo et al. 2019), only
appears at $|Z|>2\,\mathrm{kpc}$. The bending motion derived by outer-disk
stars is aligned with the recent findings that an upward bending motion exists
at $R>R_{0}$ (e.g. Wang et al. 2018; Carrillo et al. 2019; López-Corredoira et
al. 2020; Gaia Collaboration et al. 2021); by comparison, the downward bending
motion inside $R_{0}$ is trivial. The directions for both the breathing and
bending motions are mostly invariant with respect to $Z$, despite
$Z$-dependent amplitudes.

The amplitude of the breathing motion depends on both $Z$ and $R$. The spiral
density waves could drive alternate expanding and compressive motions along
the Galactocentric radius near the plane (Debattista 2014; Faure et al. 2014;
Ghosh et al. 2020), and contribute to an increasing amplitude of breathing
motions with distance to the plane within the range $|Z|\lesssim
1\,\mathrm{kpc}$ (Ghosh et al., 2020). However, considering that the observed
$|V_{\rm{breath}}|$ increases with the vertical height towards $|Z|\sim
3\,\mathrm{kpc}$, it is unlikely that the breathing mode motion is totally
attributed to spirals. As for the bending motion, Ghosh et al. (2020) argued
that the spirals are probably not responsible for the bending mode. Khoperskov
et al. (2019) proposed that the bending-mode oscillation could be induced by a
buckling bar inhabited in the disk. The differences in vertical bulk
velocities between different $\rm{[Fe/H]}$ are trivial, indicating that
different populations respond similarly to the perturbations. In this case, we
infer that both the breathing and bending motions are produced or partly
produced by long-lived external perturbations such as a satellite galaxy. It
has been argued that the perturbed vertical motions could be generated by the
Sagittarius dwarf-Milky Way interaction (Gómez et al. 2013; Widrow et al.
2014; Antoja et al. 2017; Laporte et al. 2019; Wang et al. 2019).

The north-south asymmetries in velocity dispersions are detected in this work,
though the differences between $Z>0$ stars and $Z<0$ stars are less pronounced
than those for the mean velocities. It is interesting to find that the
$\sigma_{Z}/\sigma_{R}$ in the south disk is larger than that in the north.
The velocity dispersions are probes for the secular heating process in the
disk, which are principally generated by spirals and giant molecular clouds
(GMCs). The variation in $\sigma_{Z}/\sigma_{R}$ is a measure of the variation
of the importance of the GMC as heating agents relative to spirals, since
spirals heat mainly within the plane while GMCs heat both in-plane and
vertically (Jenkins 1992; Aumer et al. 2016; Mackereth et al. 2019). The
increasing $\sigma_{Z}/\sigma_{R}$ with increasing $|Z|$ indicates that the
relative strength of spiral to GMC heating decreases with the distance to the
plane within $|Z|\lesssim 1.5\,\mathrm{kpc}$. Meanwhile, the larger
$\sigma_{Z}/\sigma_{R}$ below the plane than above suggests either a weaker
spiral heating agent or a larger number of GMCs at $Z<0$. Moreover, the
interactions between a satellite galaxy and the disk could also heat the disk
and cause north-south differences in velocity dispersions (Mackereth et al.,
2019).

The tilt of the velocity ellipsoid measured at $R<9\,\mathrm{kpc}$ is
basically consistent with the classical relation case that the
$\alpha_{\rm{tilt}}$ is nearly antisymmetric about the plane, which was
uncovered with various populations in the solar neighborhood (e.g. Binney
2014; Büdenbender et al. 2015; Ding et al. 2019; Everall et al. 2019). The
$\alpha_{\rm{tilt}}$ measured between $R=5$ and 9 kpc confirms the findings of
Hagen et al. (2019) that the orientation of the velocity ellipsoid is more
aligned with spherical when we move closer to the GC. In order to compare the
velocity ellipsoid orientation above and below the plane, we estimate the
gradients of $\alpha_{\rm{tilt}}$ with respect to $\arctan(Z/R)$ at $Z>0$ and
$Z<0$ separately. We found a gradient of $0.803\pm 0.026$, $0.784\pm 0.019$
and $0.718\pm 0.033$ at $5<R<7\,\mathrm{kpc}$, $7<R<8\,\mathrm{kpc}$ and
$8<R<9\,\mathrm{kpc}$ respectively for the north disk, as well as $0.852\pm
0.078$, $0.802\pm 0.066$ and $0.321\pm 0.075$ respectively for the south,
suggesting that the flattening effect of the tilt along the Galactocentric
radius is stronger in the south disk than in the north. The difference in the
orientation between north and south disks is more pronounced in the Galactic
anticenter. At $R>9\,\mathrm{kpc}$, when we move from $|Z|\lesssim
1\,\mathrm{kpc}$ to higher north disk, the orientation of the velocity
ellipsoid stays between a horizontal one and a spherical one. The reversed
tilt observed at $Z<-1\,\mathrm{kpc}$, evidenced by populations with different
$\rm{[Fe/H]}$, indicates that the shapes of gravitational potential for the
north and south parts of the outer disk are different when we move far away
from the plane.

The observed vertical patterns of velocity ellipsoid help us to revisit the
nature of the disk flaring. We find that the flaring of the disk begins near
the solar position, which agrees with the findings of Gyuk et al. (1999). The
position where the flaring starts to manifests itself is not settled. For
instance, López-Corredoira et al. (2002) found a strong flare beginning well
inside $R_{0}$ based on 2MASS red clump giants. To the contrary, López-
Corredoira & Molgó (2014) and Yu et al. (2021) found that the flare begins at
$R\gtrsim 10\,\mathrm{kpc}$ by using samples of SDSS-SEGUE F8V-G5V stars and
LAMOST OB stars respectively. In recent studies, the flaring of a disk has
been detected not only for our home galaxy but also for the external galaxies
(Borlaff 2016; Pinna et al. 2019; Kasparova et al. 2020), which means that a
flared disk could be a common structure in disk galaxies. The drivers of the
formation of disk flaring is still unclear, one hypothesis is the minor
mergers (Bournaud, 2009).

The vertical structures of disk kinematics found in this work call for a
deeper exploration of the origin of the long-lived north-south asymmetries in
the disk. Besides, the shapes of the vertical patterns of velocity moments as
functions of Galactocentric radius invite us to investigate the distribution
of stellar kinematics in more dimensions in the future work. Additional
measurements such as the data of chemical abundance provided by the ongoing
medium-resolution survey of the LAMOST will give further insights into the
Galactic structure and evolution.

##  6 Summary

In this work, we present an investigation into the vertical structure of
stellar velocity moments at a Galactocentric radius of $R=5-15\,\mathrm{kpc}$
and 3 kpc above and below the Galactic plane in fine detail, based on K giant
stars sampled from the LAMOST DR8. The distances of stars are estimated from
the spectroscopically derived stellar parameters of LAMOST using the technique
proposed by Carlin et al. (2015). The line-of-sight velocities are obtained
from the LAMOST catalog, and the proper motions from Gaia EDR3. The
metallicity of the sampled stars ranges from $\rm{[Fe/H]}=-1$ to 0.5. The main
results can be summarized as follows:

  * •

In the inner disk, stars both in the south and north disks are moving towards
the Galactic anticenter, with mean velocities $\langle V_{R}\rangle$
increasing with the increasing distance to the plane. In the outer disk, there
is an outward radial flow in the south disk and alternate outward and inward
flows in the north; the $\langle V_{R}\rangle$ decreases generally with the
increasing $R$, and the $\langle V_{R}\rangle$-$Z$ gradients become shallower
for larger $R$. In the range of $8<R<12\,\mathrm{kpc}$, metal-rich populations
have lower $\langle V_{R}\rangle$ than metal-poor ones.

  * •

There is contraction-like breathing mode and upward bending mode motions in
the outer disk. The amplitude of the breathing motion increases with the
increasing $|Z|$ as a whole, up to around $10\,\mathrm{km\,s}^{-1}$, and
decreases towards the Galactic anticenter. The amplitude of the bending motion
is typically within $5\,\mathrm{km\,s}^{-1}$, and nearly vanishes when $|Z|$
approaches 3 kpc. By comparison, there is only weak or trivial rarefaction-
like breathing mode and downward bending mode motions in the inner disk.

  * •

The mean azimuthal velocity decreases with the increasing distance to the
plane. When we move from the inner disk to the outer disk, the gradient
$|d\langle V_{\phi}\rangle/dZ|$ decreases from $\sim
27-28\,\mathrm{km\,s^{-1}\,kpc^{-1}}$ to $\sim
5\,\mathrm{km\,s^{-1}\,kpc^{-1}}$. In the range of $5<R<12\,\mathrm{kpc}$,
stars below the plane are rotating faster than stars above the plane. The
$\langle V_{\phi}\rangle$-$\rm{[Fe/H]}$ trend is positive at
$5<R<7\,\mathrm{kpc}$, and turning negative at $13<R<15\,\mathrm{kpc}$.

  * •

The velocity dispersions increase with the increasing distance to the plane
and decrease with the increasing $R$. The $\sigma_{Z}/\sigma_{R}$-$Z$ trend is
positive within $|Z|\lesssim 1.5\,\mathrm{kpc}$ and hardly dependent on $R$.
The $\sigma_{Z}/\sigma_{R}$ is on average larger at $Z<0$ than at $Z>0$. The
tilt of the velocity ellipsoid is more aligned to a spherical orientation at
smaller $R$, the gradient of $\alpha_{\rm{tilt}}$ with respect to
$\arctan(Z/R)$ increasing from $\sim 0.52$ at $8<R<9\,\mathrm{kpc}$ to $\sim
0.83$ at $5<R<7\,\mathrm{kpc}$. Within $9<R<15\,\mathrm{kpc}$, the velocity
ellipsoid orientation lies somewhere between a horizontal one and a spherical
one for $-1<Z<3\,\mathrm{kpc}$, while the vertical behavior of
$\alpha_{\rm{tilt}}$ reverses at $Z\lesssim-1\,\mathrm{kpc}$.

  * •

The scaleheight of the disk increases from $\lesssim 0.5\,\mathrm{kpc}$ at
$R<8\,\mathrm{kpc}$ to $\gtrsim 1.5\,\mathrm{kpc}$ at $R>12\,\mathrm{kpc}$,
indicating a flaring feature on both sides of the disk. The south disk have
larger $h_{z}$ than the north disk within $5<R<12\,\mathrm{kpc}$.

This study is supported by the National Natural Science Foundation of China
under grants Nos. 11988101, 11873052, 11890694, 11773033, 11833004 and
National Key R&D Program of China No. 2019YFA0405500. The LAMOST FELLOWSHIP is
supported by Special Funding for Advanced Users, budgeted and administrated by
Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS).
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. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences. This work
has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## Appendix

The mean and dispersion of $V_{R}$ and $V_{Z}$ are estimated based on the two-
dimensional Gaussian distribution. Let $\bm{V}=(V_{R},\,V_{Z})^{T}$ to be the
observed velocity, and $\langle\bm{V}\rangle=(\langle V_{R}\rangle,\,\langle
V_{Z}\rangle)^{T}$ to be the mean velocity. The probability density function
(PDF) of $\bm{V}$ for one star is

$f(\bm{V})=(2\pi)^{-1}\left|\bm{C_{V}}\right|^{-\frac{1}{2}}\exp\left[-\frac{1}{2}(\bm{V}-\langle\bm{V}\rangle)^{T}\bm{C_{V}}^{-1}(\bm{V}-\langle\bm{V}\rangle)\right].$
(4)

The covariance $\bm{C_{V}}$ is the sum of the dispersion tensor $\bm{D_{V}}$
and the covariance from the propagated observational uncertainties:

$\displaystyle\bm{D_{V}}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\sigma^{2}_{R}&\sigma_{RZ}\\\
\sigma_{RZ}&\sigma^{2}_{Z}\end{array}\right),$ (7) $\displaystyle\bm{C_{V}}$
$\displaystyle=$ $\displaystyle\bm{D_{V}}+\left(\begin{array}[]{cc}\delta
V_{R}^{2}&\rho\delta V_{R}\delta V_{Z}\\\ \rho\delta V_{R}\delta V_{Z}&\delta
V_{Z}^{2}\end{array}\right),$ (10)

where $\delta V_{R}$, $\delta V_{Z}$, and $\rho$ are the uncertainties in
$V_{R}$, $V_{Z}$, and their correlation coefficient. For a group of $N$ stars,
the log-likelihood distribution function for $\bm{V}$ is

$L(\langle\bm{V}\rangle,\,\sigma_{R},\,\sigma_{Z},\,\sigma_{RZ})=\sum\limits_{i=1}^{N}\ln
f_{i}(\bm{V}).$ (11)

The mean and dispersion of $\bm{V}$ can be estimated using the maximum-
likelihood method.

Figure 23: The fits to the distribution of $V_{\varphi}$ for stars at
$9<R<10\,\mathrm{kpc}$ and $1<Z<2\,\mathrm{kpc}$. The black, red, orange,
olive, and blue dots/lines denote the metallicity range $-1<\rm{[Fe/H]}<-0.6$,
$-0.6<\rm{[Fe/H]}<-0.4$, $-0.4<\rm{[Fe/H]}<-0.2$, $-0.2<\rm{[Fe/H]}<0$,
$0<\rm{[Fe/H]}<0.5$, respectively.

Unlike $V_{R}$ and $V_{Z}$, the density function of $V_{\phi}$ is non-
Gaussian. Binney (2014) proposed an analytic PDF for $V_{\phi}$, which has
been used in Ding et al. (2019) to estimate velocity moments for dwarf stars.
The PDF of $V_{\phi}$, $f(V_{\varphi})$, is written as

$\displaystyle f(V_{\varphi})$ $\displaystyle=$
$\displaystyle\rm{constant}\times\exp\big{(}-\frac{(V_{\varphi}-b_{0})^{2}}{2\sigma^{\prime
2}_{\varphi}}\big{)},$ $\displaystyle\sigma^{\prime}_{\varphi}$
$\displaystyle=$ $\displaystyle b_{1}+b_{2}V_{\varphi 100}+b_{3}V^{2}_{\varphi
100}+b_{4}V^{3}_{\varphi 100},$ (12)

where $V_{\varphi 100}\equiv V_{\varphi}/100\,\mathrm{km\,s}^{-1}$. Since the
normalizing constant is not available (Binney, 2014), we estimate
$b_{0}-b_{4}$ coupled with the normalizing constant by fitting the PDF to the
observed $V_{\varphi}$ distribution. Figure 23 gives an example of the fits
for populations with different metallicities. The mean and dispersion are
derived by taking the first and second moments of the PDF.

## References

  * Aihara et al. (2011) Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29. doi:10.1088/0067-0049/193/2/29
  * Antoja et al. (2017) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360. doi:10.1038/s41586-018-0510-7
  * Astraatmadja & Bailer-Jones (2016) Astraatmadja, T. L. & Bailer-Jones, C. A. L. 2016, ApJ, 832, 137. doi:10.3847/0004-637X/832/2/137
  * Aumer et al. (2016) Aumer, M., Binney, J., & Schönrich, R. 2016, MNRAS, 462, 1697. doi:10.1093/mnras/stw1639
  * Bailer-Jones et al. (2015) Bailer-Jones, C. A. L. 2015, PASP, 127, 994. doi:10.1086/683116
  * Bennett & Bovy (2019) Bennett, M. & Bovy, J. 2019, MNRAS, 482, 1417. doi:10.1093/mnras/sty2813
  * Binney (2014) Binney, J., Burnett, B., Kordopatis, G., et al. 2014, MNRAS, 439, 1231. doi:10.1093/mnras/stt2367
  * Borlaff (2016) Borlaff, A., Eliche-Moral, M. C., Beckman, J., et al. 2016, A&A, 591, L7. doi:10.1051/0004-6361/201628868
  * Bournaud (2009) Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJ, 707, L1. doi:10.1088/0004-637X/707/1/L1
  * Büdenbender et al. (2015) Büdenbender, A., van de Ven, G., & Watkins, L. L. 2015, MNRAS, 452, 956. doi:10.1093/mnras/stv1314
  * Carlin et al. (2013) Carlin, J. L., DeLaunay, J., Newberg, H. J., et al. 2013, ApJ, 777, L5. doi:10.1088/2041-8205/777/1/L5
  * Carlin et al. (2015) Carlin, J. L., Liu, C., Newberg, H. J., et al. 2015, AJ, 150, 4. doi:10.1088/0004-6256/150/1/4
  * Carrillo et al. (2018) Carrillo, I., Minchev, I., Kordopatis, G., et al. 2018, MNRAS, 475, 2679. doi:10.1093/mnras/stx3342
  * Carrillo et al. (2019) Carrillo, I., Minchev, I., Steinmetz, M., et al. 2019, MNRAS, 490, 797. doi:10.1093/mnras/stz2343
  * Cropper et al. (2019) Cropper, M., Katz, D., Sartoretti, P., et al. 2018, A&A, 616, A5. doi:10.1051/0004-6361/201832763
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197. doi:10.1088/1674-4527/12/9/003
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Debattista (2014) Debattista, V. P. 2014, MNRAS, 443, L1. doi:10.1093/mnrasl/slu069
  * Ding et al. (2019) Ding, P.-J., Zhu, Z., & Liu, J.-C. 2019, AJ, 158, 247. doi:10.3847/1538-3881/ab5111
  * Everall et al. (2019) Everall, A., Evans, N. W., Belokurov, V., et al. 2019, MNRAS, 489, 910. doi:10.1093/mnras/stz2217
  * Faure et al. (2014) Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564. doi:10.1093/mnras/stu428
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Antoja, T., McMillan, P., et al. 2021, arXiv:2101.05811
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2. doi:10.1051/0004-6361/201629512
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv:2012.01533
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 616, A11. doi:10.1051/0004-6361/201832865
  * Gao et al. (2015) Gao, H., Zhang, H.-W., Xiang, M.-S., et al. 2015, Research in Astronomy and Astrophysics, 15, 2204. doi:10.1088/1674-4527/15/12/008
  * Ghosh et al. (2020) Ghosh, S., Debattista, V. P., & Khachaturyants, T. 2020, arXiv:2009.02343
  * Gómez et al. (2013) Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159. doi:10.1093/mnras/sts327
  * Gravity Collaboration et al. (2018) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15. doi:10.1051/0004-6361/201833718
  * Guiglion et al. (2015) Guiglion, G., Recio-Blanco, A., de Laverny, P., et al. 2015, A&A, 583, A91. doi:10.1051/0004-6361/201525883
  * Gyuk et al. (1999) Gyuk, G., Flynn, C., & Evans, N. W. 1999, ApJ, 521, 190. doi:10.1086/307542
  * Hagen et al. (2019) Hagen, J. H. J., Helmi, A., de Zeeuw, P. T., et al. 2019, A&A, 629, A70. doi:10.1051/0004-6361/201935264
  * Jenkins (1992) Jenkins, A. 1992, MNRAS, 257, 620. doi:10.1093/mnras/257.4.620
  * Jing et al. (2016) Jing, Y., Du, C., Gu, J., et al. 2016, MNRAS, 463, 3390. doi:10.1093/mnras/stw2230
  * Karaali et al. (2014) Karaali, S., Bilir, S., Ak, S., et al. 2014, PASA, 31, e013. doi:10.1017/pasa.2014.5
  * Kasparova et al. (2020) Kasparova, A. V., Katkov, I. Y., & Chilingarian, I. V. 2020, MNRAS, 493, 5464. doi:10.1093/mnras/staa611
  * Katz et al. (2019) Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205. doi:10.1051/0004-6361/201833273
  * Khoperskov et al. (2019) Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6. doi:10.1051/0004-6361/201834707
  * Kunder et al. (2017) Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75. doi:10.3847/1538-3881/153/2/75
  * Laporte et al. (2019) Laporte, C. F. P., Minchev, I., Johnston, K. V., et al. 2019, MNRAS, 485, 3134. doi:10.1093/mnras/stz583
  * Liu et al. (2014) Liu, C., Deng, L.-C., Carlin, J. L., et al. 2014, ApJ, 790, 110. doi:10.1088/0004-637X/790/2/110
  * Lindegren et al. (2020) Lindegren, L., Bastian, U., Biermann, M., et al. 2020, arXiv:2012.01742
  * Lindegren et al. (2020) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2020, arXiv:2012.03380
  * López-Corredoira et al. (2014) López-Corredoira, M. 2014, A&A, 563, A128. doi:10.1051/0004-6361/201423505
  * López-Corredoira & Molgó (2014) López-Corredoira, M. & Molgó, J. 2014, A&A, 567, A106. doi:10.1051/0004-6361/201423706
  * López-Corredoira et al. (2002) López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., et al. 2002, A&A, 394, 883. doi:10.1051/0004-6361:20021175
  * López-Corredoira et al. (2020) López-Corredoira, M., Garzón, F., Wang, H.-F., et al. 2020, A&A, 634, A66. doi:10.1051/0004-6361/201936711
  * Luo et al. (2015) Luo, A.-L., Zhao, Y.-H., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1095. doi:10.1088/1674-4527/15/8/002
  * Mackereth et al. (2019) Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176. doi:10.1093/mnras/stz1521
  * Michalik et al. (2015) Michalik, D., Lindegren, L., & Hobbs, D. 2015, A&A, 574, A115. doi:10.1051/0004-6361/201425310
  * Monari et al. (2014) Monari, G., Helmi, A., Antoja, T., et al. 2014, A&A, 569, A69. doi:10.1051/0004-6361/201423666
  * Monari et al. (2016) Monari, G., Famaey, B., & Siebert, A. 2016, MNRAS, 457, 2569. doi:10.1093/mnras/stw171
  * Moni Bidin et al. (2012) Moni Bidin, C., Carraro, G., Méndez, R. A., et al. 2012, ApJ, 751, 30. doi:10.1088/0004-637X/751/1/30
  * Pinna et al. (2019) Pinna, F., Falcón-Barroso, J., Martig, M., et al. 2019, A&A, 625, A95. doi:10.1051/0004-6361/201935154
  * Reid & Brunthaler (2020) Reid, M. J. & Brunthaler, A. 2020, ApJ, 892, 39. doi:10.3847/1538-4357/ab76cd
  * Recio-Blanco et al. (2014) Recio-Blanco, A., de Laverny, P., Kordopatis, G., et al. 2014, A&A, 567, A5. doi:10.1051/0004-6361/201322944
  * Roeser et al. (2010) Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440. doi:10.1088/0004-6256/139/6/2440
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525. doi:10.1086/305772
  * Schönrich & Aumer (2017) Schönrich, R. & Aumer, M. 2017, MNRAS, 472, 3979. doi:10.1093/mnras/stx2189
  * Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829. doi:10.1111/j.1365-2966.2010.16253.x
  * Siebert et al. (2011) Siebert, A., Williams, M. E. K., Siviero, A., et al. 2011, AJ, 141, 187. doi:10.1088/0004-6256/141/6/187
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163. doi:10.1086/498708
  * Smith et al. (2012) Smith, M. C., Whiteoak, S. H., & Evans, N. W. 2012, ApJ, 746, 181. doi:10.1088/0004-637X/746/2/181
  * Soubiran et al. (2003) Soubiran, C., Bienaymé, O., & Siebert, A. 2003, A&A, 398, 141. doi:10.1051/0004-6361:20021615
  * Steinmetz et al. (2006) Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645. doi:10.1086/506564
  * Tian et al. (2015) Tian, H.-J., Liu, C., Carlin, J. L., et al. 2015, ApJ, 809, 145. doi:10.1088/0004-637X/809/2/145
  * Wang et al. (2019) Wang, C., Huang, Y., Yuan, H.-B., et al. 2019, ApJ, 877, L7. doi:10.3847/2041-8213/ab1fdd
  * Wang et al. (2019) Wang, H.-F., Carlin, J. L., Huang, Y., et al. 2019, ApJ, 884, 135. doi:10.3847/1538-4357/ab4204
  * Wang et al. (2018) Wang, H., López-Corredoira, M., Carlin, J. L., et al. 2018, MNRAS, 477, 2858. doi:10.1093/mnras/sty739
  * Wang et al. (2020) Wang, H.-F., López-Corredoira, M., Huang, Y., et al. 2020, MNRAS, 491, 2104. doi:10.1093/mnras/stz3113
  * Weinberg et al. (1991) Weinberg, M. D. 1991, ApJ, 373, 391. doi:10.1086/170059
  * Widrow et al. (2014) Widrow, L. M., Barber, J., Chequers, M. H., et al. 2014, MNRAS, 440, 1971. doi:10.1093/mnras/stu396
  * Widrow et al. (2012) Widrow, L. M., Gardner, S., Yanny, B., et al. 2012, ApJ, 750, L41. doi:10.1088/2041-8205/750/2/L41
  * Williams et al. (2013) Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101. doi:10.1093/mnras/stt1522
  * Wojno et al. (2018) Wojno, J., Kordopatis, G., Steinmetz, M., et al. 2018, MNRAS, 477, 5612. doi:10.1093/mnras/sty1016
  * Yanny et al. (2009) Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377. doi:10.1088/0004-6256/137/5/4377
  * York et al. (2000) York, D. G., Adelman, J., Anderson, J. E., et al. 2000, AJ, 120, 1579. doi:10.1086/301513
  * Yuan et al. (2013) Yuan, H. B., Liu, X. W., & Xiang, M. S. 2013, MNRAS, 430, 2188. doi:10.1093/mnras/stt039
  * Yu et al. (2021) Yu, Y., Wang, H.-F., Cui, W.-Y., et al. 2021, arXiv:2102.00731
  * Zhao et al. (2006) Zhao, G., Chen, Y.-Q., Shi, J.-R., et al. 2006, Chin. J. Astron. Astrophys., 6, 265. doi:10.1088/1009-9271/6/3/01
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 723. doi:10.1088/1674-4527/12/7/002