# Quantifying chemical and kinematical properties of Galactic disks

Guozhen Hu (胡国真) Key Laboratory for Research in Galaxies and Cosmology,
Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, China  University of Chinese Academy of Sciences, 19A
Yuquan Road, 100049, Beijing, China  Zhengyi Shao (邵正义) Key Laboratory for
Research in Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese
Academy of Sciences, 80 Nandan Road, Shanghai 200030, China  Key Lab for
Astrophysics, Shanghai 200234, China

###### Abstract

We aim to quantify the chemical and kinematical properties of the Galactic
disks with a sample of 119,558 giant stars having abundances and 3D velocities
taken or derived from the APOGEE DR17 and Gaia EDR3 catalogs. The Gaussian
Mixture Model is employed to distinguish the high-$\alpha$ and low-$\alpha$
sequences along the metallicity by simutaneously using the chemical and
kinematical data. Four disk components are identified and quantified that
named as h$\alpha$mp, h$\alpha$mr, l$\alpha$mp, and l$\alpha$mr disks, which
correspond to the features of high-$\alpha$ or low-$\alpha$, and metal-poor or
metal-rich. Combined with the spatial and stellar age information, we confirm
that they are well interpreted in the two-infall formation model. The first
infall of turbulent gas quickly forms the hot and thick h$\alpha$mp disk with
consequent thinner h$\alpha$mr and l$\alpha$mr disks. Then the second gas
accretion forms a thinner and outermost l$\alpha$mp disk. We find that the
inside-out and upside-down scenario does not only satisfy the overall Galactic
disk formation of these two major episodes, but also presents in the formation
sequence of three inner disks. Importantly, we reveal the inverse Age-[M/H]
trend of the l$\alpha$mr disk, which means its younger stars are more metal-
poor, indicating that the rejuvenate gas from the second accretion gradually
dominates the later star formation. Meanwhile, the recently formed stars
convergence to [M/H]$\sim$-0.1 dex, demonstrating a sufficiently mixture of
gas from two infalls.

Galaxy: disk — Galaxy: evolution — Galaxy: kinematics and dynamics — Galaxy:
structure — methods: data analysis — stars: abundances

††software: Astropy (Astropy Collaboration et al., 2013), Numpy (van der Walt
et al., 2011), Scipy (Oliphant, 2007), Matplotlib (Hunter, 2007), emcee
(Foreman-Mackey et al., 2013).

##  1 Introduction

The Milky Way (MW) is a unique laboratory to study galaxy formation and
evolution, where stellar properties can be measured and analyzed in great
detail. Galactic archaeology utilizes stellar abundances and kinematics as
relics to investigate the history of the MW, which relies on the assumption
that any history of the Galactic formation and dynamical evolution will leave
its marks on chemical abundances and kinematics, as well as their
correlations, of stellar populations (Freeman & Bland-Hawthorn, 2002).

In investigating the assembly history of the MW, it is essential to separate
different components (i.e., stellar populations). The traditional view holds
that there are two disk components, a thin disk and a thick disk in the MW.
Gilmore & Reid (1983) revealed this pattern about 40 years ago, and their work
showed that at the Galactic South pole, the vertical stellar density profile
could not be described by only one exponential. Then, the traditional
definition of thin and thick disk features have been widely used: the thick
disk is older than the thin disk (e.g., Wyse & Gilmore 1988; Haywood et al.
2013; Bensby et al. 2014; Bashi & Zucker 2019); has a larger rotational lag
compared to the thin disk (e.g., Guiglion et al. 2015; Rojas-Arriagada et al.
2016; Belokurov et al. 2020; Bashi et al. 2020); and has a higher velocity
dispersion than the thin disk (e.g.,Recio-Blanco et al. 2014; Guiglion et al.
2015; Wojno et al. 2016).

Earlier studies used kinematics to separate thick and thin disk populations in
two different ways. The first is to assume that the space velocities of
different populations have different three-dimensional Gaussian distributions
(Soubiran & Girard, 2005; Bensby & Feltzing, 2006; Ruchti et al., 2010; Katz
et al., 2011; Kordopatis et al., 2011). The second is to distinguish different
components directly based on the Tommre diagram (Re Fiorentin et al., 2019).
However, as mentioned by previous studies (e.g., Reddy et al. 2006; Bensby et
al. 2014), those distinctions based on kinematical criteria are imperfect.
That is because the thin and thick disks overlap in their kinematical
distributions, and the heating or migration processes might enhance such a
mixture. Then, it might misidentify transitional stars, which perhaps bias
conclusions (Bovy et al., 2012).

With the accumulation of high-resolution spectroscopic data from large survey
projects, studies show that the $\alpha$-enhancement ([$\alpha$/M]) and the
metallicity ([M/H]) of local stars follow two distinct tracks (e.g., Adibekyan
et al. 2011; Bensby et al. 2014), which provide much information about stellar
populations, such as the star formation efficiency, the accretion history, and
the gas outflows (e.g., Andrews et al. 2017; Weinberg et al. 2017). The
$\alpha$-elements are produced mainly during the Type II supernovae (SNe)
explosions on a short time scale, whereas other elements (e.g., the iron) are
mainly produced by Type Ia SNe on a much longer time scale. So the
[$\alpha$/M] can be used as a ’clock’ to probe the star formation history of a
Galactic component (Schuster et al., 2012). The high-$\alpha$ population is
typically associated with the thick disk, while the low-$\alpha$ stars are
associated with the thin disk. Moreover, many works realize that the
dissecting method of components based on chemistry is superior to the method
based on kinematics, since the chemical abundance is a relatively stable
property of stars, which will not evolve after their formation. Thus, the
mono-abundance populations (MAPs) are commonly used to trace the MW disk
formation processes.

At the metal-poor part ([M/H] $<-0.2$ dex), the stars are clearly separated in
terms of [$\alpha$/M]. However, in the metal-rich range, the distinction is
less obvious (Fuhrmann, 2004; Adibekyan et al., 2011). Barbuy et al. (2018)
also revealed that the high-$\alpha$ and low-$\alpha$ sequences are merged up
to [M/H] $\sim$ 0.15 dex. Some works suggested that the thick disk might not
include stars with [M/H] $>-0.3$ dex (Reddy et al., 2006). Usually, metal-rich
stars with relatively higher [$\alpha$/M] are regarded as a particular
population. Sometimes, the metal-rich range of the high-$\alpha$ sequence is
called the ’transition zone’ (Duong et al. 2018, see their Fig. 8), or is
classified as the so called h$\alpha$mr disk (e.g. Adibekyan et al. 2011,
2013, see their Fig. 1; Yu et al. 2021, see their Fig. 3), or is named as a
’bridge’ component (Ciucă et al., 2021). In the metal-rich range, most
previous works artificially distinguish high- or low-$\alpha$ components
(thick or thin disks) with straight lines on [$\alpha$/M]-[M/H] plane (e.g.,
Lee et al. 2011; Anders et al. 2014; Re Fiorentin et al. 2019). Definitely,
such a simple cut may bias the properties of the real h$\alpha$mr components
and also contaminate the thin disk at the metal-rich part. Therefore, it is
valuable to accurately determine the disks in the metal-rich region, if there
are surely two separate populations.

The velocity and chemical abundance of stars are thought to involve clues to
the formation and evolution of the MW. Stellar populations with different
kinematics and abundances are related to different formation epochs. The
origin of thick disk is still a matter of debate. A possible scenario is that
the thick disk stars are formed by the dynamical heating of a pre-existing
thin disk through small satellite mergers (e.g., Kazantzidis et al. 2008;
Schönrich & Binney 2009a; Sharma et al. 2021) or/and directly accreted from
disrupted satellites (relatively metal-poor) around the MW (Abadi et al. 2003;
Bird et al. 2013). Another scenario is based on the two-infall model (e.g.,
Chiappini et al. 1997; Grisoni et al. 2017; Spitoni et al. 2019; Lian et al.
2020; Spitoni et al. 2020). In this model, the thick disk is due to an earlier
gas infall which occurred a strong star burst and a subsequent secular star
forming process. Then, a second infall rejuvenates the metal-poor gas to form
the kinematically colder and thin disk. Anyway, for both of these scenarios,
more detailed constraints should be based on the quantitative description of
disks with their kinematical and chemical properties, and their coincidence.
Therefore, it is necessary to investigate the separation and connection of
disk components by simultaneously using kinematical and chemical data.

Some works have already combined the kinematical and chemical information in
separating the thin and thick disks. Duong et al. (2018) dissected disks by
constructing a mixture model in the space of [$\alpha$/M], [M/H] and radial
velocity. They obscured the proper motion due to the noteworthy observational
errors at that time. In Nikakhtar et al. (2021), they fitted a mixture of a
number of Gaussian distributions using the three dimensional velocities and
iron abundances of stars, while leaving the number of components as a free
parameter. Unfortunately, these works have not adopted the full information of
kinematics or chemical abundance.

In this work, we propose to quantitatively describe the disk chemical and
kinematical properties with a sample of giant stars from the Apache Point
Observatory Galactic Evolution Experiment (APOGEE) survey. Combining with the
astrometric data from the Gaia, the sample stars have chemical measurements
([M/H] and [$\alpha$/M]) and three-dimensional velocity data. Then, for
different [M/H] bins, we employ the Gaussian Mixture Model (GMM) to fit the
star distribution in the remained four-dimensional data space, i.e., one for
the [$\alpha$/M] and the other three for the velocities. This approach will
allow us to curve the stellar population properties across the whole
metallicity range and rigorously solve the mean values and their dispersions
of these separated stellar components. It is worthwhile to mention that there
is a potential advantage of using velocity and abundance distributions, since
both of them are hardly affected by the incompleteness of the survey area, or
the observational selection effects.

Additionally, based on the GMM framework, we can calculate the membership
probabilities for individual stars of each component, i.e., belonging to the
high-$\alpha$ or the low-$\alpha$ sequences. These membership probabilities
are also valuable in further statistical analysis with other variables, e.g.,
the stellar age or their Galactocentric distances.

The paper is organized as follows. In Section 2, we describe the giant star
sample used in this work, taken from the APOGEE and Gaia catalogs, supply with
the derived distance, velocities, and stellar age measurements. Section 3
explains the strategy and method that we used to fit the data and assess the
results. In Section 4, we present the fitting results leading to the
classification of four disk components, the membership of individual stars,
and discuss the effectiveness of our method. Then, We quantitatively describe
the chemical abundance, kinematics, spatial and age distributions for each
disk component and discuss the implications of the disk formation scenario in
Section 5. Conclusions are summarized in Section 6. Throughout this work, we
use [M/H] as the measure of metallicity and [$\alpha$/M] as a tracer of
$\alpha$-enhancement.

##  2 DATA AND SAMPLE SELECTION

Our sample is assembled by cross-matching APOGEE DR17 (Abdurro’uf et al.,
2021) and Gaia Early Data Release 3 (EDR3) (Gaia Collaboration et al., 2021;
Lindegren et al., 2021). The sample was then restricted to giant stars to
avoid potential differential trends in the chemical measurements of dwarfs and
giants. From the APOGEE data, we take the [M/H], [$\alpha$/M], and radial
velocity of stars. The [M/H] and [$\alpha$/M] values are derived from the
ASPCAP analysis of combined spectrum, and the $\alpha$-elements considered are
O, Mg, Si, S, Ca and Ti. Combining with the coordination, Gaia proper motion,
and distance (Bailer-Jones et al., 2021) derived from the Gaia parallax, we
calculate each star’s position and velocity coordinates (the 6D phase space)
and subsequently transfer them to the Galactocentric coordinate system.
Moreover, we also collect stellar age measurements from different sources for
further discussions.

###  2.1 APOGEE DR17

APOGEE DR17 is a part of the final data release of the SDSS-IV (Sloan Digital
Sky Survey; Blanton et al. 2017) via APOGEE-1 and APOGEE-2 programs. It
contains 733,901 AllStar entries of 656,940 individual stars across both North
and South hemispheres (Abdurro’uf et al., 2021).

Firstly, in order to avoid the bias of stellar element abundances obtained by
ASPCAP when dealing with different types of stars (Ahumada et al., 2020), in
this work, we only consider the giant stars (ASPCAP_GRID = g). Secondly, since
we are only interested in the general structure of the MW and mainly focus on
the disks, we exclude stars with PROGRAMNAME = bulge, gc, oc and halo-dSph.
After that, the data quality is further controlled, and we only consider the
stars with ASPCAP_CHI2 $<10$, ASPCAP_FLAG $=0$, SNR $>70$, and M_H_ERR $<0.1$
dex. Thus, we have totally picked up 235,250 giants from APOGEE DR17. For the
duplicate observed sources, we use the error-weighted means instead.

###  2.2 Gaia EDR3

Gaia EDR3 contains $\sim$1.81 billion sources with the magnitude completeness
to $G\sim 20.7$ mag and has parallaxes and proper motions for a subset of
$\sim$ 1.47 billion sources (Gaia Collaboration et al., 2021; Lindegren et
al., 2021). Gaia EDR3 extended the observational baseline to 34 months
resulting in more precise astrometric measurements. Compared with the Gaia
DR2, the parallax error of EDR3 was reduced by a factor of $\sim$20$\%$, and
the accuracy of the proper motion was increased twice.

We cross-match the sample of giants from APOGEE DR17 with Gaia EDR3. We remain
the stars with epsi $<4$ and sepsi $<8$ (Soltis et al., 2021) to ensure the
sources with acceptable astrometric solutions. In order to desert astrometric
binaries and other anomalous sources, we retain only those stars with the re-
normalized unit weight error RUWE $<$ 1.4 (Lindegren et al., 2021).

Finally, 119,558 stars are left in our sample. Fig. 1 shows the distribution
of sample stars in the chemical abundance plane of [$\alpha$/M] versus [M/H].
Generally, this diagram is applied to recognize the chemically defined thick
(high-$\alpha$) and thin (low-$\alpha$) disk populations. The bimodal sequence
of [$\alpha$/M] is very significantly shown by our sample. More interestingly,
there is a low number density area close to [M/H]$\simeq-0.1$ and
[$\alpha$/M]$\simeq 0.14$, which implies a break of the star formation
process.

###  2.3 Distance and velocities

The coordination and proper motion of sample stars are directly taken from
Gaia EDR3. We employ the photo-geometric distance estimated by Bailer-Jones et
al. (2021) to avoid the parallax bias and high uncertainties of the Gaia
targets far away from the Sun. Combined with the APOGEE’s radial velocity, one
can easily calculate the 6D kinematic information of a single star (3D
positions and 3D velocities) related to the Sun. Then, we transform them to
the Galactocentric cylindrical system $(R,\phi,Z,V_{R},V_{\phi},V_{Z})$ by
assuming that the Galactocentric distance of the Sun is $R_{\odot}=8.125$ kpc
(Gravity Collaboration et al., 2018) and its height from the Galactic plane is
$Z_{\odot}=20.8$ pc (Bennett & Bovy, 2019). We also use a solar motion of
(11.1, 242, 7.25) km s-1 in the radial, rotational and vertical directions,
respectively (Schönrich et al. 2010; Bennett & Bovy 2019).

Figure 1: Number density distribution of sample stars in the
[$\alpha$/M]-[M/H] plane. The colored distribution is for the full sample, and
the dashed-lines correspond to the 1, 2 and 3$\sigma$ contours of stars with
Age-S18 value.

###  2.4 Age measurements

The accurate stellar age determinations are decisive to investigate the MW
formation and evolution (Marín-Franch et al., 2009; Soderblom, 2010; Minchev
et al., 2019). Unfortunately, stellar age is challenging to be measured
directly (e.g., Wyse & Gilmore 1988; Haywood et al. 2013; Bensby et al. 2014).
It can only be estimated through model-dependent (e.g., Jørgensen & Lindegren
2005; Takeda et al. 2007) or empirical methods (e.g., Barnes, 2007).

In this work, we employ two recently published age catalogs. The first one we
considered is from Sanders & Das (2018) (hereafter Age-S18). Based on the
Bayesian framework, they used the photometric, spectroscopic, and astrometric
information to characterize the stellar parameters. This catalog provides age
estimation for $\sim$ 3 million stars. We reject stars with age uncertainties
$\sigma_{\tau}>$ 1.5 Gyr. Then, 86137 APOGEE/Gaia/Age-S18 stars are left. In
Fig. 1, the dashed contour lines draw the distribution of the sample with
Age-S18, which coincides with the colored map, indicating that the Age-S18
subset has no obvious selection bias.

The second age we used was taken from the SDSS DR17 Value-Added Catalog
(Mackereth et al., 2019) (hereafter Age-M19). Age-M19 is estimated with a
Bayesian neural network model (Leung & Bovy, 2019) trained by the
asteroseismic ages (Pinsonneault et al., 2018). Mackereth et al. (2019)
claimed that the Age-M19 of stars with [Fe/H] $<$ -0.5 dex are less reliable
than the others because of the lack of training set stars at the metal-poor
region. The typical uncertainty of this age is 0.1 dex (25 percent). Totally,
there are 118,766 APOGEE/Gaia/Age-M19 stars remaining in our sample.

##  3 Methodology

To model the distribution of our sample stars, we first divide them in [M/H]
bins. Then, for each bin, we use the multivariate GMM to fit the subsample
with four variables, [$\alpha$/M], $V_{R},V_{\phi}$ and $V_{Z}$. This binning
strategy is according to the following factors. Firstly, as Fig. 1
demonstrated, the [M/H] distributions are not regular (both high- and
low-$\alpha$ sequences show ’banana’ shapes), so they cannot be directly
modeled as simple profiles. Secondly, [M/H] covers a wide range, and compared
with [$\alpha$/M], has much smaller observational errors. Moreover, in the
most part of the [$\alpha$/M]-[M/H] plane, for a given [M/H] value, the
[$\alpha$/M] appears clearly bimodal shape. We suppose that each
$\alpha$-enhancement sequence corresponds to a specific star formation
process. Thus, for a given [M/H], the distribution of [$\alpha$/M] of each
sequence is well expected to have a Gaussian profile (equal to the log-normal
distribution of $\alpha/{\rm M}$), due to its complex origin through the star
forming physical conditions. Therefore, it will be very convenient to analyze
the distribution structure with such binned subsamples.

###  3.1 Gaussian Mixture Model

A GMM is a parametric probability density function represented as a weighted
sum of Gaussian components for multivariate densities. For the specific
subsample within a given bin of [M/H], we assume that the stars follow a
weighted sum of a number ($K$) of four dimensional (4D) normal distributions.
If we further assume that this 4D normal distribution is independent of our
adopting variables, [$\alpha$/M] and three velocities
($V_{R},V_{\phi},V_{Z}$), then it could be fully shaped with the parameters of
mean values and dispersions for each one-dimensional normal distribution in
each variable space, i.e.,
$({\boldsymbol{\mu}},{\boldsymbol{\sigma}})=(\mu_{[\alpha/{\rm
M]}},\sigma_{[\alpha/{\rm
M}]};\mu_{V_{R}},\sigma_{V_{R}};\mu_{V_{\phi}},\sigma_{V_{\phi}};\mu_{V_{Z}},\sigma_{V_{Z}})$.

For the $i$th star with observational data $\boldsymbol{x}_{i}=([\alpha/{\rm
M}]_{i},V_{R,i},,V_{\phi,i},V_{Z,i})$, the likelihood of this star belonging
to the GMM is written as:

$\mathcal{L}_{i}=\sum_{k=1}^{K}w_{k}\boldsymbol{\mathcal{N}}_{k}(\boldsymbol{x}_{i}|\boldsymbol{\mu}_{k},\boldsymbol{\sigma}_{k}),$
(1)

where the weights of each components satisfy

$\sum_{k=1}^{K}w_{k}=1.$ (2)

Then we can write down the joint logistical likelihood for the subsample with
$N$ data points,

$\log\mathcal{L}=\sum_{i=1}^{N}\log\mathcal{L}_{i}.$ (3)

We employ the EM algorithm (Isobe et al., 1986) method to search the
$\mathcal{L}_{\rm max}$ across all parameter ranges to obtain their optimal
values.

It is clear that we need to somehow decide on the number ($K$) of Gaussian
components by trading off quality of fit against the number of introduced free
parameters. To do so, we use the Bayesian Information Criterion (BIC; Schwarz
1978; Connolly et al. 2000) to determine how many mixtures we should use. The
BIC is defined as:

${\rm BIC}=-2\log\mathcal{L}_{\rm max}+m\log(N),$ (4)

where $m$ is the number of free parameters, and in this work, $m=3K-1$, for
each sub-samples. Then, we compare the corresponding BICs with different
numbers of components and select the model with the smallest BIC value.

###  3.2 Membership probability and effectiveness

For a set of optimal parameters of
$(w_{k},\boldsymbol{\mu}_{j},\boldsymbol{\sigma}_{k};k=1\,{\rm to}\,K)$, the
probability that a data point ($\boldsymbol{x}_{i}$ for the $i$th star)
belonging to the $j$th population is:

$P_{j}(i)=\frac{w_{j}\boldsymbol{\mathcal{N}}_{j}(\boldsymbol{x}_{i}|\boldsymbol{\mu}_{j},\boldsymbol{\sigma}_{j})}{\sum_{k=1}^{K}w_{k}\boldsymbol{\mathcal{N}}_{k}(\boldsymbol{x}_{i}|\boldsymbol{\mu}_{k},\boldsymbol{\sigma}_{k})}.\\\
$ (5)

Alternatively, one can calculate the $P_{j}$ value in the case of involving
only a part of observational variables. For example, if we neglect the star’s
velocities and consider the [$\alpha$/M] data only, we can determine the GMM
distributions in the [$\alpha$/M] space and calculate the corresponding
probability $P_{j}^{\alpha}$. Similarly, one can also calculate probability
$P_{j}^{\boldsymbol{v}}$ as the probability that only using the velocity data.
These different definitions of $P_{j}$ will be used to evaluate the
improvement of the combination of multiple observational data (see Sec. 4.4
for some examples).

It is clear that $P_{j}$ is a kind of normalized values through $K$
components. Based on the distribution of $P_{j}$ for the $j$th component, we
can quantitatively evaluate how effective the results of our membership
determination are (Shao & Zhao, 1996). The effectiveness of membership
determination is defined as:

$E_{j}=1-\frac{N\sum_{i=1}^{N}{P_{j}(i)[1-P_{j}(i)]}}{\sum_{i=1}^{N}P_{j}(i)\sum_{i=1}^{N}[1-P_{j}(i)]}.$
(6)

The range of $E$ value is $[0,1]$. The bigger $E$ is, the more effective the
membership determination is. In other words, $E$ could be regarded as an index
to assess the mixture level of the GMM, e.g., a more discrete Gaussian
component will lead to a bigger $E$ value from its membership probability
distribution.

##  4 Identification of Stellar Components and their fitting results

###  4.1 Number of mixed components

All of our sample stars have been binned in 0.1 dex interval of [M/H].
Totally, there are 28 bins from -2.3 to 0.5 dex of [M/H]. In this work, we
need to choose one or two components to model the data within a given bin. So,
we employ either one or two components GMM to fit each subsample and calculate
the corresponding values of BIC1 or BIC2, and make the choice of $K$ based on
$\Delta{\rm BIC}=({\rm BIC}_{1}-{\rm BIC}_{2})/N$. The $\Delta$BIC results are
listed in table 1. Among them, there are 16 bins having $\Delta{\rm BIC}>0$,
all within the range from [M/H]=$-1.2$ to 0.4 dex. So, for these bins, the two
Gaussian distributions ($K=2$) are adopted. According to the mean values of
[$\alpha$/M], they are named as high- or low-$\alpha$ populations. For the
other bins, at the metal-poor region ([M/H]$<-1.2$ dex) and the most metal-
rich bin ([M/H]$\simeq 0.45$ dex), their $\Delta{\rm BIC}<0$, which suggests
it is unnecessary to separate the subsample into two populations.

It should be mentioned that, at the metal-rich part ([M/H]$>-0.1$ dex),
although the bimodal feature of [$\alpha$/M] is not as clear as those at the
more metal-poor region, the BIC values indicate that only one Gaussian is not
enough to curve the distributions for most of these bins. Obviously, this is
an evidence of none-single population from the data themselves. It implies the
complicated composition of these metal-rich stars. In other words, the so-
called h$\alpha$mr disk might exist, but definitely could not be simply
distinguished by straight lines on the [$\alpha$/M]-[M/H] plane.

###  4.2 Classification and parameters of disk components

Figure 2: Results of the Gaussian mixture decomposition. The top-left panel is
the projection on the [$\alpha$/M]-[M/H] plane. Other three panels are
projections on three velocity-[M/H] planes respectively. The symbols with
error bars show the GMM fitting results of the mean and dispersion values for
each Gaussian component, with ’red’ and ’black’ for the bi-normal model
($K=2$), while ’green’ for the single Gaussian model ($K=1$). Different
symbols represent different populations, with diamonds for the halo, squares
and circles for the high- and low-$\alpha$ sequences respectively, and the
only green star for the most metal-rich sub-sample. Contours show the
probability-weighted number densities of high-$\alpha$ (red) and low-$\alpha$
(black) stars. The solid lines correspond the 1, 2 and 3$\sigma$ contours, and
the black dashed line is the 4$\sigma$ contour of the low-$\alpha$ sequence to
demonstrate the different distribution features of halo and disk stars. The
blue dashed vertical lines at [M/H]=-0.1 divide the metal-poor and metal-rich
parts.

The fitting results of parameters $(\boldsymbol{\mu},\boldsymbol{\sigma})$ are
summarized in table 1. The mean and dispersion values are also shown as the
symbols with error bars in Fig. 2. Since there are higher fitting precision,
i.e., the uncertainties of all parameters are significantly smaller than the
corresponding Gaussian dispersions, we omit the errors of fitting parameters
either in the table or in the figure.

Based on these chemical and velocity parameters, we can classify each Gaussian
component of the MW populations. For the metal-poor subsamples with
[M/H]$<-1.2$ dex, $K=1$ for each bins. They have relatively higher
[$\alpha$/M], roughly zero means and large dispersions ($\sim 100$ km s-1) of
all velocity components, $V_{R}$, $V_{\phi}$ and $V_{Z}$. Obviously, they are
halo stars. In the range of $-1.2<$[M/H]$<0.4$, where we should use two
Gaussian components to model the subsamples, it can be further discriminated
into several parts. For the high-$\alpha$ sequence, it could be split at
[M/H]$\sim-0.1$ dex. The metal-poor part of this sequence is regarded as the
canonical thick disk, and we use h$\alpha$mp to denote this component. For the
metal-rich part of this sequence, the [$\alpha$/M] further decreases with
metallicity, but still higher than another component in the same [M/H] bin. As
noted by previous works (Adibekyan et al. 2011, 2013; Yu et al. 2021), we
define this component as the h$\alpha$mr disk.

The low-$\alpha$ sequence also has multiple compositions. The first split
point is at [M/H]$\simeq-0.8$, where there is a clear break of the kinematic
features, either in the mean value of $V_{\phi}$, or in all three velocity
dispersions. The more metal-poor bins follow the kinematics of the halo
component, though their [$\alpha$/M] values are lower than the corresponding
thick disk stars. This phenomenon could be explained as the MW halo is a
complicated population that contains quite a lot of accreted dwarf galaxies
(Nissen & Schuster, 2010). The other more metal-rich bins of this sequence are
similar, with their characteristics representing the canonical thin disk,
either in chemical abundance or in kinematics. However, here we prefer to
separate these bins into two groups by [M/H]$\sim-0.1$. The component with
$-0.8<$[M/H]$<-0.1$ of the low-$\alpha$ sequence is regarded as the
l$\alpha$mp disk, while the more metal-rich component ([M/H]$>-0.1$) is named
as the l$\alpha$mr disk. According to the fitting parameters, there are slight
differences between these two low-$\alpha$ disks. The l$\alpha$mp disk shows a
correction between [$\alpha$/M] and [M/H], with the more metal-rich stars
having lower [$\alpha$/M] values. Whereas the [$\alpha$/M] values of the
l$\alpha$mr disk do not change with the metallicity. On the other hand, the
$\sigma_{V_{Z}}$ of the l$\alpha$mr disk is smaller than the l$\alpha$mp disk.
Rigorously speaking, there are overlaps between two low-$\alpha$ disks.
However, considering the break in the [$\alpha$/M]-[M/H] correlation and other
evidences from the age and spatial distributions (see detailed comparison and
discussions in Sec. 5 ), we use a simple cut at [M/H]=-0.1 to separate these
two populations.

For the most metal-rich bin of our sample ([M/H]$\sim$0.45), the optimal
distribution model is a single Gaussian. We suppose it is a merge of the
h$\alpha$mr and l$\alpha$mr populations.

Table 1: The fitting results and mean ages in different metallicity bins

[M/H] bin | N |  $\Delta$BIC | K | $f_{h\alpha}$ |  Comp. a |  [$\alpha$/M] | $\sigma_{[\alpha/M]}$ | $V_{R}$ | $\sigma_{V_{R}}$ | $V_{\phi}$ | $\sigma_{V_{\phi}}$ | $V_{Z}$ | $\sigma_{V_{Z}}$ |  Age-S18 b |  Age-M19 c
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
dex | - | - | - | - | - | dex | dex | ${\rm km\,s^{-1}}$ | ${\rm km\,s^{-1}}$ | ${\rm km\,s^{-1}}$ | ${\rm km\,s^{-1}}$ | ${\rm km\,s^{-1}}$ | ${\rm km\,s^{-1}}$ | Gyr | Gyr
[-2.3, -2.2] | 15 | -0.73 | 1 | … | halo | 0.323 | 0.050 | -7.2 | 70.6 | -2.0 | 98.7 | 11.4 | 109.7 | 10.85 | 8.85
[-2.2, -2.1] | 41 | -0.30 | 1 | … | halo | 0.295 | 0.081 | 11.6 | 106.5 | -6.7 | 110.0 | 24.9 | 107.5 | 11.03 | 8.74
[-2.1, -2.0] | 63 | -0.27 | 1 | … | halo | 0.311 | 0.082 | 20.0 | 102.3 | 27.9 | 90.7 | -1.8 | 113.7 | 11.10 | 8.86
[-2.0, -1.9] | 89 | -0.32 | 1 | … | halo | 0.289 | 0.081 | -2.1 | 127.5 | 50.1 | 93.0 | -0.3 | 97.8 | 10.80 | 8.83
[-1.9, -1.8] | 116 | -0.33 | 1 | … | halo | 0.269 | 0.081 | 16.9 | 122.9 | 24.4 | 99.5 | -5.2 | 96.4 | 10.84 | 8.53
[-1.8, -1.7] | 140 | -0.24 | 1 | … | halo | 0.247 | 0.068 | 3.4 | 136.6 | 53.0 | 98.8 | 0.4 | 92.4 | 10.80 | 8.46
[-1.7, -1.6] | 179 | -0.08 | 1 | … | halo | 0.221 | 0.080 | 36.5 | 126.6 | 65.4 | 98.3 | -10.7 | 110.6 | 10.28 | 7.98
[-1.6, -1.5] | 263 | -0.12 | 1 | … | halo | 0.227 | 0.074 | -1.4 | 124.4 | 48.6 | 92.2 | -18.9 | 98.3 | 10.40 | 7.84
[-1.5, -1.4] | 294 | -0.13 | 1 | … | halo | 0.228 | 0.062 | 2.6 | 127.4 | 41.2 | 84.3 | -25.0 | 97.0 | 9.05 | 7.51
[-1.4, -1.3] | 359 | -0.10 | 1 | … | halo | 0.226 | 0.065 | -13.8 | 137.1 | 39.0 | 84.1 | -20.3 | 99.7 | 8.86 | 6.98
[-1.3, -1.2] | 441 | -0.01 | 1 | … | halo | 0.226 | 0.064 | -19.8 | 142.4 | 38.8 | 76.8 | -2.0 | 107.8 | 8.40 | 6.24
[-1.2, -1.1] | 418 | 0.08 | 2 | 0.32 |  h$\alpha$mp | 0.289 | 0.062 | 2.4 | 82.0 | 133.9 | 76.4 | 3.4 | 70.8 | 8.83 | 7.20
|  |  |  |  | halo | 0.210 | 0.092 | -4.8 | 126.8 | 24.6 | 80.4 | -10.2 | 95.2 | 8.90 | 5.82
[-1.1, -1.0] | 422 | 0.02 | 2 | 0.38 |  h$\alpha$mp | 0.284 | 0.056 | 0.2 | 81.8 | 140.7 | 73.2 | -1.3 | 66.3 | 8.63 | 6.90
|  |  |  |  | halo | 0.185 | 0.095 | -5.2 | 122.9 | 30.4 | 82.5 | -13.0 | 94.8 | 8.85 | 5.96
[-1.0, -0.9] | 502 | 0.18 | 2 | 0.52 |  h$\alpha$mp | 0.284 | 0.059 | -3.5 | 81.9 | 142.9 | 70.2 | 0.9 | 60.3 | 8.49 | 7.68
|  |  |  |  | halo | 0.143 | 0.041 | 6.3 | 128.8 | 31.0 | 87.4 | -2.3 | 92.8 | 8.07 | 5.37
[-0.9, -0.8] | 883 | 0.57 | 2 | 0.74 |  h$\alpha$mp | 0.282 | 0.074 | 4.6 | 80.3 | 142.9 | 73.0 | 1.4 | 59.0 | 8.17 | 7.51
|  |  |  |  | halo | 0.130 | 0.063 | 2.7 | 114.4 | 85.9 | 106.4 | 7.0 | 84.2 | 6.66 | 4.37
[-0.8, -0.7] | 2040 | 0.64 | 2 | 0.86 |  h$\alpha$mp | 0.263 | 0.056 | -0.5 | 79.0 | 136.4 | 76.1 | -0.1 | 54.1 | 8.10 | 7.81
|  |  |  |  |  l$\alpha$mp | 0.098 | 0.025 | -4.4 | 28.0 | 216.2 | 19.6 | 2.2 | 28.0 | 5.93 | 4.82
[-0.7, -0.6] | 4223 | 1.04 | 2 | 0.79 |  h$\alpha$mp | 0.262 | 0.052 | -1.4 | 73.9 | 142.9 | 71.7 | 1.7 | 51.5 | 7.79 | 8.16
|  |  |  |  |  l$\alpha$mp | 0.097 | 0.023 | 0.6 | 32.3 | 218.7 | 21.2 | 3.1 | 26.2 | 5.33 | 5.07
[-0.6, -0.5] | 8485 | 1.59 | 2 | 0.50 |  h$\alpha$mp | 0.258 | 0.045 | -3.2 | 64.7 | 170.6 | 55.5 | -2.0 | 44.7 | 7.46 | 8.30
|  |  |  |  |  l$\alpha$mp | 0.101 | 0.025 | 5.7 | 36.4 | 221.7 | 21.6 | 2.5 | 29.8 | 5.30 | 5.01
[-0.5, -0.4] | 14110 | 1.62 | 2 | 0.36 |  h$\alpha$mp | 0.246 | 0.038 | 9.3 | 65.7 | 179.3 | 45.7 | -3.6 | 44.4 | 7.07 | 8.42
|  |  |  |  |  l$\alpha$mp | 0.093 | 0.026 | 0.6 | 34.8 | 225.8 | 19.9 | -1.5 | 27.3 | 4.83 | 4.47
[-0.4, -0.3] | 15871 | 1.48 | 2 | 0.31 |  h$\alpha$mp | 0.232 | 0.034 | 0.7 | 58.9 | 180.9 | 45.0 | -0.9 | 39.1 | 6.89 | 8.46
|  |  |  |  |  h$\alpha$mp | 0.078 | 0.026 | 3.2 | 34.0 | 220.6 | 22.3 | 2.4 | 20.2 | 4.39 | 3.78
[-0.3, -0.2] | 16683 | 1.20 | 2 | 0.27 |  h$\alpha$mp | 0.205 | 0.040 | 0.7 | 52.3 | 187.3 | 41.2 | -1.7 | 37.3 | 6.86 | 8.21
|  |  |  |  |  l$\alpha$mp | 0.060 | 0.026 | -0.0 | 35.1 | 219.6 | 21.6 | 1.2 | 18.1 | 4.03 | 3.40
[-0.2, -0.1] | 16336 | 0.91 | 2 | 0.25 |  h$\alpha$mp | 0.140 | 0.056 | 2.4 | 51.3 | 196.6 | 35.0 | 1.2 | 30.5 | 6.37 | 7.44
|  |  |  |  |  l$\alpha$mp | 0.040 | 0.023 | -1.0 | 33.1 | 218.3 | 20.6 | 1.0 | 15.2 | 3.54 | 3.08
[-0.1,  0.0] | 14137 | 0.59 | 2 | 0.25 |  h$\alpha$mr | 0.098 | 0.050 | -3.2 | 47.5 | 199.7 | 30.5 | 1.8 | 27.8 | 6.14 | 7.12
|  |  |  |  |  l$\alpha$mr | 0.026 | 0.022 | -0.7 | 33.0 | 217.8 | 21.5 | -0.4 | 14.9 | 3.57 | 3.41
[0.0,  0.1] | 11257 | 0.36 | 2 | 0.27 |  h$\alpha$mr | 0.071 | 0.040 | -1.4 | 47.3 | 204.7 | 30.7 | 0.8 | 26.2 | 5.82 | 6.64
|  |  |  |  |  l$\alpha$mr | 0.016 | 0.022 | -3.3 | 33.7 | 216.3 | 21.6 | 0.3 | 14.8 | 3.59 | 4.08
[0.1,  0.2] | 7270 | 0.21 | 2 | 0.18 |  h$\alpha$mr | 0.062 | 0.031 | 2.2 | 42.6 | 202.4 | 30.3 | 3.0 | 28.1 | 5.63 | 6.60
|  |  |  |  |  l$\alpha$mr | 0.013 | 0.019 | 0.3 | 35.6 | 214.2 | 23.1 | -1.1 | 15.5 | 4.08 | 4.65
[0.2,  0.3] | 3490 | 0.15 | 2 | 0.17 |  h$\alpha$mr | 0.050 | 0.033 | -1.4 | 47.2 | 199.4 | 32.6 | -3.6 | 32.2 | 5.44 | 6.33
|  |  |  |  |  l$\alpha$mr | 0.016 | 0.023 | 0.4 | 36.7 | 211.6 | 23.9 | -1.3 | 16.0 | 4.57 | 5.35
[0.3,  0.4] | 1310 | 0.38 | 2 | 0.08 |  h$\alpha$mr | 0.019 | 0.038 | -2.2 | 46.9 | 201.3 | 32.0 | -0.5 | 25.1 | 5.07 | 5.78
|  |  |  |  |  l$\alpha$mr | 0.017 | 0.024 | 1.4 | 34.9 | 209.2 | 23.9 | 1.2 | 16.0 | 5.08 | 5.78
[0.4,  0.5] | 121 | -0.38 | 1 | … | mr | 0.015 | 0.020 | 4.5 | 36.7 | 203.0 | 34.8 | -1.4 | 16.5 | 4.81 | 5.90

  * a

Classification of the MW components.

  * b

Sanders & Das (2018). The Age-S18 is probability-weighted average age of each
component in each bin.

  * c

Mackereth et al. (2019). The Age-M19 is probability-weighted average age of
each component in each bin. The Age-M19 of stars with [M/H]$<-0.5$ dex are
less reliable.

###  4.3 Membership probability of each component

Using Eq. 5, for the bimodal bins, we calculate the probability of each star
belonging to the high-$\alpha$ sequence, $P_{h\alpha}$. Accordingly, the
low-$\alpha$ probability, $P_{l\alpha}=1-P_{h\alpha}$. The probability values
are listed in table 2. Based on the disk classification, we can obtain the
probabilities of each stars belonging to a specific disk or halo components.
By summing of the probabilities, in our giant star sample, the fraction of
different components are, halo (2.5%), h$\alpha$mp disk (24.7%), h$\alpha$mr
disk (7.3%), l$\alpha$mp disk (41.4%) and l$\alpha$mr disk (24.1%).

Table 2: Chemical and kinamatical data and membership probability of sample
stars

APOGEE-ID a | SOURCE-ID b | [$\alpha$/M] c | [M/H] d | $V_{R}$ | $V_{\phi}$ | $V_{Z}$ | $K$ | $P_{h\alpha}$
---|---|---|---|---|---|---|---|---
|  | dex | dex | km s-1 | km s-1 | km s-1 |  |
2M10061549+4720367 | 822188376907896064 | 0.245 | -2.247 | -127.6 | 113.1 | -64.8 | 1 | …
2M15505257+2756175 | 1223770877101987712 | 0.279 | -1.855 | 203.9 | 2.4 | 181.7 | 1 | …
2M21024981-0118487 | 6917139903604771712 | 0.286 | -0.929 | -64.2 | -2.1 | -77.0 | 2 | 0.99
2M23472814+6149358 | 2012869825245864320 | 0.123 | -0.607 | 0.4 | 200.1 | -45.6 | 2 | 0.02
2M02325239+2829337 | 128147809635150336 | 0.134 | -0.462 | -54.2 | 147.8 | -8.7 | 2 | 0.90
2M07032260+5715085 | 1000544349938145920 | 0.109 | -0.205 | -98.6 | 198.0 | -24.0 | 2 | 0.35
2M18461632-2259448 | 4078443430416791040 | 0.077 | -0.199 | -29.4 | 243.9 | 53.0 | 2 | 0.82
2M14512752-5925456 | 5878295071153144192 | 0.041 | 0.140 | -42.5 | 225.9 | -9.8 | 2 | 0.10
2M07060606-7440158 | 5261235195577098880 | 0.018 | 0.255 | 45.8 | 214.3 | -40.7 | 2 | 0.30
2M15181135-4920358 | 5902178284791783808 | 0.008 | 0.491 | 16.7 | 176.9 | 0.2 | 1 | …

  * a

APOGEE DR17 object name.

  * b

Gaia EDR3 source id.

  * c

[$\alpha$/M] in APOGEE DR17.

  * d

[M/H] in APOGEE DR17.

Figure 2 also shows the probability weighted number densities of high- and
low-$\alpha$ sequences in red or black contours. As expected, in the
[$\alpha$/M]-[M/H] plane, these two sequences distinguish clearly in the
region of more metal-poor than -0.1 dex, but still with unneglectable overlap.
The overlapping phenomena are more significant in those $V$-[M/H] planes than
in the [$\alpha$/M]-[M/H] plane. However, it shows differences in the velocity
dispersions and the mean values of $V_{\phi}$. These kinematical differences
enhance the distinctness of these two [$\alpha$/M] dominated sequences. In
all, these two sequences have natural looking density profiles in either
[$\alpha$/M]-[M/H] or $V$-[M/H] planes, which implies the combination of
chemical and kinematic data is valuable.

Additionally, the probability values provide an approach to take into account
all individual stars’ contribution to a specific population. This illustrates
the ability of mixture modeling to allow full probabilistic analysis, which is
particularly significant when considering the degree of overlap between
populations.

###  4.4 Effectiveness of the combination of abundance and kinematics

In Fig. 3, as examples, we plot the fitting results and number density
distributions in the [$\alpha$/M]-$V_{\phi}$ plane of three representational
subsamples ([M/H]$\simeq-0.95$, -0.35 and 0.15 dex, for cases of the
h$\alpha$mp-halo, h$\alpha$mp-l$\alpha$mp and h$\alpha$mr-l$\alpha$mr
mixtures), to demonstrate the reasonableness of using the bimodal distribution
in the [$\alpha$/M] versus velocity space. Besides the most significantly
discrete h$\alpha$mp-l$\alpha$mp distribution (in the middle panel), the other
two subsamples are also shown indications of the bimodality in the
[$\alpha$/M]-$V_{\phi}$ plane.

The importance of the combination of chemical and kinematical data can also be
shown in Fig. 4, in which we plot the probability distributions of these three
subsamples. For each subsample, there are three kinds of probabilities, with
$P_{h\alpha}^{\alpha}$ for only using the [$\alpha$/M] data,
$P_{h\alpha}^{\boldsymbol{v}}$ for only using the velocities, and
$P_{h\alpha}^{\alpha\boldsymbol{v}}$ for simultaneously using both of them.
Respectively, the effectiveness indexes $E$ are calculated through Eq. 6, and
denoted as $E^{\alpha}$, $E^{\boldsymbol{v}}$ and $E^{\alpha\boldsymbol{v}}$.
Comparing with these $E$ values, we can find that, the [$\alpha$/M] is the
most effective variable to separate the h$\alpha$mp-l$\alpha$mp mixture (in
the middle panel with the largest value of $E^{\alpha}$), while the kinematic
data are very effective to separate the halo and h$\alpha$mp components (the
left panel). In all, when both [$\alpha$/M] and velocity data are involved,
the effectiveness of the separation are increased, more or less.

We can also calculate and compare the $E$ values for the whole high- and
low-$\alpha$ sequences. The corresponding values are: $E^{\alpha}=0.88$,
$E^{\boldsymbol{v}}=0.67$ and $E^{\alpha\boldsymbol{v}}=0.93$. It shows the
improvement on the combination of chemical and kinematical data, and also
demonstrate the high effectiveness of the membership determination of the
current sample.

Figure 3: The number density and GMM fitting results for three
representational metallicity bins ([-1, -0.9], [-0.4, -0.3], [0.1, 0.2]) in
the $V_{\phi}$\- [$\alpha$/M] plane. Red and black symbols with error bars
show the mean and dispersion values of each Gaussian component. Red and Black
lines are 1 and 2$\sigma$ contours of the probability-weighted number density
for each component.  Figure 4: Probability distributions of the high-$\alpha$
components of three representational metallicity bins (same as in Fig. 3). The
green, blue and red histograms represent the membership probabilities using
only velocity or chemical data, or both of them. The effectiveness indexes for
each probability distribution are calculated.

##  5 Comparison and Discussion

Based on the chemical and kinematical data and the GMM method, the MW disk is
divided into four components with distinctive chemical properties,
high-$\alpha$ or low-$\alpha$, metal-poor or metal-rich, named as h$\alpha$mp,
h$\alpha$mr, l$\alpha$mp and l$\alpha$mr disks respectively. Here we describe
their properties individually, discuss the correlations among them and further
analyze the implications for the formation scenario of the MW disk. The MW
halo is also a separated component characterized by its chemical abundance and
kinematics, but will not be discussed in this paper.

###  5.1 Chemical Abundance

h$\alpha$mp disk

It is the major part of the high-$\alpha$ sequence, covering the [M/H] range
from -1.2 to -0.1 dex, with monotonically decreased [$\alpha$/M] from 0.29 to
0.14 dex. This significant [$\alpha$/M]-[M/H] anti-correlation could be
reproduced by the chemical evolution model of a star forming process. Along
this [$\alpha$/M]-[M/H] trend, the [$\alpha$/M] dispersion
$\sigma_{[\alpha/{\rm M}]}$ is about 0.05 dex, varies from 0.035 to 0.074 dex.

h$\alpha$mr disk

It extends the high-$\alpha$ sequence towards the most metal-rich point
[M/H]$\sim 0.4$ dex, with [$\alpha$/M] decreasing from 0.10 to $\sim$ 0.02
dex. In the transition area, the h$\alpha$mp and h$\alpha$mr disks seem to
have breaks between them for the average values of number density and their
distributions in the [$\alpha$/M]-[M/H] plane. For the main part of this disk,
$\sigma_{[\alpha/{\rm M}]}\lesssim 0.04$ dex, which is slightly smaller than
the h$\alpha$mp disk.

l$\alpha$mp disk

It covers the [M/H] from -0.8 to -0.1 dex. This population distinct from the
h$\alpha$mp disk with a clear offset of [$\alpha$/M] $\sim$ -0.15 dex, while
has the similar trend of the [$\alpha$/M]-[M/H] correlation, with [$\alpha$/M]
decreases from 0.10 to 0.04 dex. The $\sigma_{[\alpha/{\rm M}]}$ is generally
smaller than the high-$\alpha$ sequence, with the values around 0.025 dex.
However, the number density contours in the upper-left panel of Fig. 2 show a
definite overlap of these two metal-poor disks. That is because they are not
only determined by chemical abundance, but also by involving star’s
kinematical characteristics. Therefore, we declare that the dispersion of
these two populations should be more reliable.

l$\alpha$mr disk

Unlike the other three disk components, the l$\alpha$mr disk appears almost
constant mean values of [$\alpha$/M]$\sim 0.016$ dex across the [M/H] range
from -0.1 to 0.4 dex, with similar dispersions $\sigma_{[\alpha/{\rm M]}}\sim
0.020$ dex, which is the smallest value among the four components. Although
the number densities of l$\alpha$mr disk may have strong overlaps with the
h$\alpha$mr disk, it is possible to distinguish them via the mean [$\alpha$/M]
values, dispersions, and the [$\alpha$/M]-[M/H] trends.

In brief, these four disk components show different distribution features in
the [$\alpha$/M]-[M/H] plane. The dispersions of the [$\alpha$/M] are vary
among them, but generally agree with other works (e.g. $\sigma_{[\alpha/{\rm
Fe}]}\sim$ 0.04 dex in Vincenzo et al. 2021). Moreover, since the
observational error of [$\alpha$/M] is not too large, with a typical value
$\sim 0.01$ dex for stars with [M/M] $>$ -1.0 dex, one can suppose that the
[$\alpha$/M] dispersions are mainly due to the intrinsic distributions. The
largest [$\alpha$/M] dispersion of the h$\alpha$mp disk implies that it has
been under the most violent star formation process.

###  5.2 Kinematics

The common kinematic features of four disk components are their average values
($\sim 0$ km s-1) of radial velocity $V_{R}$, vertical velocity $V_{Z}$, and
the symmetrical distributions of these two velocities. Obviously, these disks
distinguish from each others through the rotational velocity $V_{\phi}$ and
the velocity dispersions of all three velocities
($\sigma_{V_{R}},\sigma_{V_{\phi}},\sigma_{V_{Z}}$), as we itemize below.

h$\alpha$mp disk

The most significant feature of this disk is the $V_{\phi}$-[M/H] correlation,
with $V_{\phi}$ increase with [M/H] from 134 to 196 km s-1. A liner fitting
arrives a positive slope of d$V_{\phi}$/d[M/H]= 49 km s-1 dex-1. This slope
value is compatible with the earlier observational or simulative studies (Lee
et al. 2011; Adibekyan et al. 2013; Kawata et al. 2018 ). On the other hand,
the velocity dispersions, including $\sigma_{V_{R}}$, $\sigma_{V_{\phi}}$ and
$\sigma_{V_{Z}}$, all monotonously decrease with [M/H] from 82 to 51 km s-1,
76 to 35 km s-1 and 71 to 31 km s-1, respectively.

h$\alpha$mr disk

Its $V_{\phi}\approx 200$ km s-1, roughly keeps in constant within the metal-
rich range. It is a little bit slower than the typical rotational speed of the
MW. Velocity dispersions are all having small variations, with the average
values of $\sigma_{V_{R}}\sim$45 km s-1, $\sigma_{V_{\phi}}\sim$32 km s-1 and
$\sigma_{V_{Z}}\sim$27 km s-1.

l$\alpha$mp disk

This disk generally obeys the typical rotation of the MW, with
$V_{\phi}\sim$220 km s-1 across its whole [M/H] coverage from -0.8 to -0.1
dex. A very small gradient with [M/H] could be found, which means the rotation
of more metal-rich stars is a few km s-1 slower. The slope is about -7 km
s-1dex-1. It is significantly shallower than those slopes of Lee et al. (2011)
whose subsample of thin disk is G-dwarfs, and Recio-Blanco et al. (2014) whose
sample of thin disk is FGK-type stars. The velocity dispersions
$\sigma_{V_{R}}\sim 33$ km s-1 and $\sigma_{V_{\phi}}\sim 22$ km s-1. While
the $\sigma_{V_{Z}}$ shows a gradient with [M/H] significantly, decreases from
28 to 15 km s-1, which means the metal-poor part (usually consists of older
stars) of this disk has larger $V_{Z}$ dispersions.

l$\alpha$mr disk

Its $V_{\phi}$ follows the small [M/H] gradients of the l$\alpha$mp disk , and
further decreases to 210 km s-1. The velocity dispersions are approximately
remained constants with $\sigma_{V_{R}}\sim 35$ km s-1, $\sigma_{V_{\phi}}\sim
22$ km s-1, $\sigma_{V_{Z}}\sim 16$ km s-1. It is the dynamically coldest
disk.

In summary, the two metal-rich disks and the l$\alpha$mp disk follow the
rotation of the MW. The h$\alpha$mp disk is slowly rotating, while the
velocity increases with the metallicity. Based on the comparison of velocity
dispersions, we can conclude that the h$\alpha$mp is the hottest disk, while
the l$\alpha$mr is the coldest disk. The other two components could be
regarded as warm disks but still are rotation-dominated populations.

It is interesting that there are breaks between the two high-$\alpha$ disks,
in the trends of $V_{\phi}$-[M/H] and $\sigma_{V}$-[M/H] for all three
velocity dispersions. For these four kinematical parameters, the values
increase/decrease with [M/H] for the h$\alpha$mp disk, but keep in constants
for the h$\alpha$mr disk. These kinds of breaks imply the change of star
formation procedure. In contrast, the kinematics of the l$\alpha$mp and
l$\alpha$mr disks are quite similar, and there are no significant breaks
between them, except the $\sigma_{V_{Z}}$ which is obviously larger of the
l$\alpha$mp disk than that of the l$\alpha$mr disk.

Additionally, one may find that, for all disks and in all [M/H] bins except
the $\sigma_{V_{Z}}$ of the most metal-poor part of the l$\alpha$mp disk, we
have $\sigma_{V_{R}}>\sigma_{V_{\phi}}>\sigma_{V_{Z}}$. It follows the general
relations of the whole MW disk from previous works (Bensby et al., 2003; Bond
et al., 2010).

Finally, it is worth noting that, for the low-$\alpha$ sequence, at the
[M/H]$\sim-0.8$ dex, there is an extremely sudden break of $V_{\phi}$ and
dispersions of all three velocities. This is a very clear separation of the
l$\alpha$mp disk and the halo, though there is a continuity in the chemical
abundance distributions.

###  5.3 Spatial distribution

Figure 5: Left panel: The cumulative $R_{\rm gc}$ distributions for the
h$\alpha$mp, h$\alpha$mr, l$\alpha$mp and l$\alpha$mr disks. Right panel: The
cumulative $|$Z$|$ distributions for the same four disks.

Our sample is limited and shaped by the APOGEE survey, and is not volume
completed. So it is difficult to accurately model the spatial distribution of
individual disk components, e.g., to measure their scale-length or scale-
height. Nevertheless, we can still compare their relative distributions. For
this purpose, we adopt the cumulative distribution function (CDF) of $R_{\rm
gc}$ and $|Z|$. Fig. 5 shows these two CDFs of the four disk components
individually. Here we accumulate the membership probability of each star that
belongs to the specific disk component.

In the left panel of Fig. 5, we can find that, the h$\alpha$mp disk is the
innermost component with its members mainly distribute from $R_{\rm gc}=$ 3.5
to 10 kpc. While the l$\alpha$mp disk stars mainly distribute farther than 7.5
kpc from the Galactic center, and extend to $\sim$15 kpc. The other two metal-
rich disks, h$\alpha$mr and l$\alpha$mr are somewhat in between of them,
whereas more close to the h$\alpha$mp disk. Therefore, we can roughly define
the h$\alpha$mp, h$\alpha$mr, and l$\alpha$mr as three inner disks, while the
l$\alpha$mp as an outer disk. Moreover, along the chemical abundance evolution
trace, i.e., the increase of [M/H] or the decrease of [$\alpha$/M], the sizes
of the three inner disks, from h$\alpha$mp to h$\alpha$mr and then to
l$\alpha$mr, enlarge in sequence.

The right panel of Fig. 5 shows the pattern of three thin disks verse one
thick disk. The l$\alpha$mr disk is the thinnest one, with its sample stars
concentrated within $\pm 1$ kpc from the Galactic plane. The thickest disk is
the h$\alpha$mp which extends its member stars to higher than 4 kpc. The other
two intermediate disks look quite similar in the $|Z|$ distributions. However,
since the l$\alpha$mp disk has a much larger $R_{\rm gc}$ distribution, it
should be even thinner than the h$\alpha$mr disk. Interestingly, considering
the three inner disks, h$\alpha$mp, h$\alpha$mr and l$\alpha$mr, they are not
only enlarged, but also getting thinner in sequence.

Combined with the $R_{\rm gc}$ and $|Z|$ distributions, we can conclude that
the h$\alpha$mp disk is the only real thick (hot) disk with the inner most
distribution towards the Galactic center. The two low-$\alpha$ disks, either
the outer or the inner ones, are thin (cold) disks. Without considering the
l$\alpha$mp disk, the other three inner disks follow the inside-out and
upside-down formation scenario (Brook et al. 2006; Bird et al. 2013; Minchev
et al. 2015). These results are also consistent with recent analysis by
Spitoni et al. (2021) using the APOGEE DR16.

###  5.4 Age distribution

Figure 6: The age-metallicity distributions. Left and right panels are for
Age-S18 (Sanders & Das, 2018) and Age-M19 (Mackereth et al., 2019)
respectively. The red and black contours represent the 1 and 2$\sigma$
probability-weighted number densities of high- and low-$\alpha$ disks,
respectively. The symbols are probability-weighted mean ages of each
metallicity bin. Green arrows in the left panel hint the evolution trace of
four disks. As mentioned by Mackereth et al. (2019), the Age-M19 of stars with
[Fe/H] $<-0.5$ dex are less reliable, hence we use hollow symbols for [M/H]
$<-0.5$ dex.

For most of our sample stars, we have the age measurements from two sources
(see Sec. 2.4 for details), Age-S18 and Age-M19. Although both of them have
large uncertainties for individual stars, it is surely possible for us to
statistically determine the typical age of a subsample within a given [M/H]
bin, and then explore the age features of four disks. It is worth noting that,
in the GMM, we only use the chemical and kinematical data, so the age could be
regarded as an independent variable in the discussion.

Fig. 6 shows the probability-weighted number densities of high- and
low-$\alpha$ sequences as contours in the Age-[M/H] plane. Age-S18 and
Age-M19, are plotted in two panels separately. The symbols show the average
age of each disk component in each [M/H] bin (also listed in table 1). Since
the Age-M19 of stars with [Fe/H] $<-0.5$ dex are less reliable than the others
(Mackereth et al., 2019), we use hollow symbols for [M/H] poorer than -0.5
dex, and neglect them in following analysis.

With the current stellar age measurements, either Age-S18 or Age-M19, large
dispersions ($\sim$1.5 Gyr) of both high- and low-$\alpha$ sequences are
present. However, there is an overall trend that metal-rich stars are
noticeably younger than metal-poor stars. It is also clear that there is a
time lag between two sequences, with the high-$\alpha$ populations roughly
2$\sim$3 Gyrs older than the low-$\alpha$ ones.

Three disk components, the h$\alpha$mp, h$\alpha$mr, and l$\alpha$mp, have
similar correlations between age and [M/H], which indicates the more metal-
rich stars are averagely younger than the more metal-poor stars. These trends
obey the normal abundance enhancement procedure of stellar formation for a
population. The age distribution of the h$\alpha$mp disk is in rough agreement
with that of the canonical thick disk stars previously discussed by many works
(e.g., Bensby et al. 2007; Anders et al. 2018; Miglio et al. 2021; Lagarde et
al. 2021).

On the contrary, and more importantly, this kind of correlation of the
l$\alpha$mr disk is inverse, with younger stars tending to be more metal-poor
than older stars. This phenomenon is clearly shown for both Age-S18 and
Age-M19. The abnormal Age-[M/H] trend might be explained as there is a
sufficient gas rejuvenation in the later forming stage of this disk. It will
cause the enhancement of the star-forming activity leading to much more stars
in the l$\alpha$mr disk than in the h$\alpha$mr disk, which is exactly the
case we can show in their fractions and number density distributions (see in
Fig. 2). It also causes an interesting phenomenon that the youngest population
of our MW is not the most metal-rich but has [M/H]$\sim$-0.1, which is a
little bit poorer than the Sun.

Alternatively, if we compare the stars of the low-$\alpha$ sequence along the
timeline, one may notice that the older stars have wide, or probably
dichotomous metallicity distribution, while the younger stars’ distribution is
gradually getting narrower and converges to [M/H]$\sim$-0.1 dex. This
phenomenon is even more obvious for the Age-M19 distribution. It implies that
in the early stage of these two low-$\alpha$ disks, their stars formed
separately in outer and inner regions by irrelevant gas with significantly
different metalicities. Then, with ongoing gas enrich and mixture processes,
the later formed stars tend to have similar metalicities.

From another point of view, the three inner disks seem to conjoint in a
consequent formation process, whereas the outer disk, l$\alpha$mp, is more
likely to occupy an independent formation process. In the outer region, it
enhances the disk itself, and also brings metal-poor gas into the inner region
to mix with metal-rich gas and then leading to the abnormal Age-[M/H]
relationship of the l$\alpha$mr disk.

###  5.5 Differences among disks

Among the four disk components, there are conjoint components, two
high-$\alpha$ disks, two low-$\alpha$ disks, and two heavily mixed metal-rich
disks. Actually, there are pieces of evidence to identify the independence of
these stellar populations, as we summarize below.

h$\alpha$mp $vs$ h$\alpha$mr

In the [$\alpha$/M]-[M/H] plane (see red contours of the top-left panel of
Fig. 2), there is a sudden decrease of [$\alpha$/M] from h$\alpha$mp to
h$\alpha$mr, and also a possible saddle point of the number density, at
[M/H]$\simeq-0.1$ and [$\alpha$/M]$\simeq 0.14$. The most significant
difference between these two disks is the change of the $V_{\phi}$ and the
break of the $[V_{\phi}]$-[M/H] correlation. Also, the velocity dispersions
$\sigma_{V}$ of the h$\alpha$mp disk are all larger than the h$\alpha$mr disk.
Another significant difference is that the $|Z|$ dispersion of h$\alpha$mp
stars are much higher than those of the h$\alpha$mr disk.

l$\alpha$mp $vs$ l$\alpha$mr

There is no break in the number density distribution of these two disks.
However, there are breaks in the [$\alpha$/M]-[M/H] and $\sigma_{V_{Z}}$-[M/H]
correlations. One of the obvious differences is in their $R_{\rm gc}$
distributions, where the l$\alpha$mp stars mostly located in the outer
regions, which may lead to the well known radial gradient of metallicity. The
most important difference is shown in the Age-[M/H] trends, where the
l$\alpha$mr disk follows an inverse abundance evolution trace with younger
stars having poorer metallicity. This is the dominant reason to distinguish
the conventional thin disk into two parts.

h$\alpha$mr $vs$ l$\alpha$mr

These two disks are directly and rigorously determined by the GMM. So the
smaller BIC value for a two components fit ($K=2$) is the first evidence of
the existence of two metal-rich disks. Besides the difference of the
[$\alpha$/M] distributions, the h$\alpha$mr disk also has lower $V_{\phi}$ and
larger dispersions for all three velocities than the l$\alpha$mr disk. Of
cause, the most important difference between them is their inverse Age-[M/H]
trends. Generally, the h$\alpha$mr stars are older than the l$\alpha$mr stars
at the same metallicity. This characteristic has been mentioned by Kordopatis
et al. (2011) and Lagarde et al. (2021).

It is commonly agreed that there is a transition zone between the thick disk
and the metal-rich part of the thin disk. This stellar family, also named as
the h$\alpha$mr disk were already reported and investigated by many previous
works with different star samples (Adibekyan et al. 2011; Ness et al. 2013;
Adibekyan et al. 2013; Kawata et al. 2018; Yu et al. 2021). It is still an
open question whether this population only represents the metal-rich tail of
the thick (h$\alpha$mp) disk, or the h$\alpha$mr stars might have originated
from the inner Galactic bulge and have migrated towards the solar
neighborhood. Unfortunately, previous works usually employed a simple cut of
higher $\alpha$-enhancement stars at the metal-rich region. While in this
work, we use the GMM to divest the h$\alpha$mr disk from the metal-rich part
of the thin (low-$\alpha$) disk, so it can rigorously derive all features of
this disk. This rigorous statistical solution also benefits the characterizing
of the remaining l$\alpha$mr component.

In summary, the h$\alpha$mr and l$\alpha$mr are two particular disks that we
have divested, identified, and quantitatively described rather than in
previous works.

###  5.6 Implication of the disk star formation scenario

Our results suggest a formation scenario of the MW disks involving two
pathways of star formation and chemical evolution, which generally supports
the two-infall picture of Milky Way (e.g. Chiappini et al. 1997; Spitoni et
al. 2020). The first pathway is from the h$\alpha$mp disk to the h$\alpha$mr
disk, and then to the l$\alpha$mr disk. It occurs the initial starburst with
its consequent process. The canonical thick (h$\alpha$mp) disk forms earlier
from the highly turbulent gas with a broad vertical spread. This episode is
very short and confirms the results of other works (Chiappini et al. 1997;
Spitoni et al. 2020). After the rapid star formation is quenched, it transfers
to a long-term secular evolution phase with lower SFR, leading to also inner
but kinematically cold disk components. In the later stage of this star-
forming episode, the remaining gas was further cooled down but continually
rejuvenated with lower [M/H] gas from the second gas accretion, resulting in
the gradually decrease of the average stellar metallicity. In the view of
kinematics, these three components satisfy a process of going from hot to cold
and accelerating rotation.

The second pathway is from the l$\alpha$mp disk to the younger part of the
l$\alpha$mr disk. It may be triggered by another accretion of a huge amount of
low [M/H] gas about 3 Gyr later than the first starburst and occurs the
formation of the outer thin disk (the l$\alpha$mp disk), and the significant
time lag confirms the results of previous works (Spitoni et al. 2019, 2020,
2021). Then, the metal-poor or somewhat enriched gas continually rejuvenated
the most metal-rich gas of the inner part, and gradually dominate the
l$\alpha$mr disk’s recent star formation. Comparing with the correlation
slopes of Age-[M/H] of the two metal-poor disks, we can find the slope of
l$\alpha$mp disk is steeper than the h$\alpha$mp disk. That means this second
star formation process may not be as violent as the earlier one. The
discrepancy of the $\alpha$-enhancement and the velocity dispersions between
these two disks also support this conclusion, since the $\alpha$-enhancement
is much more efficiency for the shorter time scale star formation.

This two pathways scenario also leads to a complicated formation process of
the l$\alpha$mr disk and reminds us that the metal-rich region ([M/H]$>$-0.1)
is far away from the detailed investigation. Although the separation of
h$\alpha$mr and l$\alpha$mr disks may not be the only deconstruction, it
provides a novel view in understanding the origin of the metal-rich stars.

This overall MW disk formation scenario reminds several interesting MAP
points. One is at [M/H]$\simeq-0.1$ and [$\alpha$/M]$\simeq 0.14$, which is
the breakpoint of two high-$\alpha$ disks. According to the age measurements
used in this paper, this phase-transfer happened at $\sim$6.5 (Age-S18) or
$\sim$7.5 (Age-M19) Gyrs ago. Another break point is between the l$\alpha$mp
disk and halo components, at [M/H]$\simeq-0.8$ and [$\alpha$/M]$\simeq 0.1$.
This point indicates the beginning of the second gas accretion happened about
6.5 Gyrs ago. Besides, the most metal-rich population of the MW disk is at
[M/H]$\simeq 0.45$ and [$\alpha$/M]$\simeq 0.0$, which is the turn-off point
showing the gas rejuvenation making effects about 5.0 (Age-S18) or 6.0
(Age-M19) Gyrs ago, and is roughly 1.5 Gyrs after the quench of the first
starburst and the beginning of the second gas infall. Moreover, the youngest
MAP of the star sample of this work is at [M/H]$\simeq-0.1$ and [$\alpha$/M]
$\simeq 0.0$, which is a little bit more metal-poor than the Sun and
represents the chemical abundance of the remained and well mixed cool gas. We
noted that, although the current age measurements for a large sample of stars
may still have large uncertainties or unneglectable bias, the relative order
and the time lag of different populations are reliable.

Our results are also consistent with the well-known inside-out and upside-down
MW disk formation picture (Brook et al. 2006; Bird et al. 2013; Minchev et al.
2015). The overall scenario is caused by the two major gas infall processes
mainly represented by two metal-poor disks. The early formed h$\alpha$mp disk
is significantly smaller and thicker (hotter) than the second accretion
occurred l$\alpha$mp (outer) disk. Besides this large scale picture, if we
only consider the three inner disks of the first episode, we can also discover
an inside-out and upside-down behavior. Along the formation sequence, from
h$\alpha$mp to h$\alpha$mr, and then to l$\alpha$mr, these inner disks are
sequentially enlarged in size and getting thinner in thickness.

We have to mention that the MW dynamical processes, e.g., the radial migration
(Schönrich & Binney, 2009a, b) or the heating process (Bournaud et al., 2009)
by the accretion of dwarf galaxies (or sub-halos), may significantly modulate
the kinematics and spatial distributions of different disk components formed
in different epochs. Therefore, a more detailed kinematical analysis,
including the information of the stellar orbit properties, energy, and angular
momentum, is required. We may extensively discuss these issues in a separate
paper.

##  6 Conclusion

This work aims to quantify the chemical and kinematical properties of the
Galactic disk components to understand their structure formation and
evolution. We use a large sample of 119,558 giants from APOGEE DR17 that
cross-match Gaia EDR3. By employing the GMM method for sub samples of
metallicity bins, we obtain the distribution parameters (mean and dispersion
values) of [$\alpha$/M], $V_{R}$, $V_{\phi}$, and $V_{Z}$ as functions of
[M/H] (table 1 and Fig. 2). So the dichotomous distribution feature is then
identified and quantified, not only in the $\alpha$-enhancement but also in
the velocity spaces. However, we still use high- and low-$\alpha$ to name the
distinct sequences which can well characterize them. As a by-product, we
calculate the membership probabilities belonging to each sequence for
individual sample stars (table 2), and the membership determination
effectiveness is also evaluated (Fig. 3 and 4).

These two sequences show significant overlaps in the distributions on
[$\alpha$/M]-[M/H] and $V$-[M/H] planes which confirm the necessity of the
combination of multiple variables in separating these two sequences.
Obviously, properties vary along the sequences, so we can determine four disk
components with distinctive chemical properties and name them as h$\alpha$mp,
h$\alpha$mr, l$\alpha$mp, and l$\alpha$mr disks, respectively. Compared with
previous works (Soubiran & Girard, 2005; Bensby & Feltzing, 2006; Adibekyan et
al., 2011; Bensby et al., 2014), the most improvement of this classification
is that we divest the h$\alpha$mr disk from the metal-rich part with a
rigorous statistical method and distinguish the l$\alpha$mr disk from the
canonical thin disk. It makes the l$\alpha$mr disk a particular and well-
defined component. Additionally, two specific MAPs are figured out as the
breakpoints between the two high-$\alpha$ disks and between the l$\alpha$mp
disk and halo populations.

Besides the difference in chemical abundance, four disk components show their
difference in the properties of kinematics, spatial distributions (Fig. 5),
$V$-[M/H] correlations (Fig. 2), or Age-[M/H] correlations (Fig. 6). The
varieties between these disks generally support the two-infall formation
scenario of the MW disk (e.g., Chiappini et al. 1997; Grisoni et al. 2017;
Spitoni et al. 2019; Lian et al. 2020). The first infall of turbulent gas
quickly forms the hot thick disk and consequently changes to the secular phase
with a lower star formation rate to form the h$\alpha$mr and l$\alpha$mr disks
(Lian et al., 2020; Spitoni et al., 2021). The second major gas accretion is
more temperate to form the l$\alpha$mp disk that dominates the outer Galactic
disk region and rejuvenate the inner disk part to influence the secular
formation process of the l$\alpha$mr disk. It results in a more metal-poor
young population.

We have discovered that the inside-out and upside-down disk forming scenario
is hierarchical. The h$\alpha$mp disk, which is the canonical thick disk that
dominates the first gas infall is significantly smaller and thicker than the
l$\alpha$mp formed with the later gas accretion (Brook et al. 2006; Bird et
al. 2013; Minchev et al. 2015). Meanwhile, the three inner disks, h$\alpha$mp,
h$\alpha$mr, and l$\alpha$mr, also show the evidence of getting extender and
thinner successively.

The most interesting phenomenon we have revealed is the inverse Age-[M/H]
trend of the l$\alpha$mr disk, which means that younger stars are more metal-
poor than older stars of this disk component. It violates the metal increasing
trend of a spatially enclosed star-forming area and indicates that, in the
inner Galactic disk region, the rejuvenate gas gradually takes the initiative
of the later star formation.

Finally, from the perspective of the entire Galactic disk, the recently formed
stars convergence to [M/H]$\sim$-0.1 dex. It demonstrates that the gas of the
second accretion is sufficiently mixed with the remains of the first gas
infall.

We thank Ruixiang Chang, Maosheng Xiang for the helpful discussions. This work
is supported by National Key R&D Program of China No. 2019YFA0405501, the
National Natural Science Foundation of China (NSFC) under grant U2031139, and
the science research grants from the China Manned Space Project with NO. CMS-
CSST-2021-A08. 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.

The stellar parameters, abundances, RVs from APOGEE DR17 were derived from the
allStar files available at https://www.sdss.org/dr17/. Funding for the Sloan
Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the
U.S. Department of Energy Office of Science, and the Participating
Institutions.

## References

  * Abadi et al. (2003) Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 591, 499, doi: 10.1086/375512
  * Abdurro’uf et al. (2021) Abdurro’uf, Accetta, K., Aerts, C., et al. 2021, arXiv e-prints, arXiv:2112.02026.  https://arxiv.org/abs/2112.02026
  * Adibekyan et al. (2011) Adibekyan, V. Z., Santos, N. C., Sousa, S. G., & Israelian, G. 2011, A&A, 535, L11, doi: 10.1051/0004-6361/201118240
  * Adibekyan et al. (2013) Adibekyan, V. Z., Figueira, P., Santos, N. C., et al. 2013, A&A, 554, A44, doi: 10.1051/0004-6361/201321520
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Anders et al. (2018) Anders, F., Chiappini, C., Santiago, B. X., et al. 2018, A&A, 619, A125, doi: 10.1051/0004-6361/201833099
  * Anders et al. (2014) —. 2014, A&A, 564, A115, doi: 10.1051/0004-6361/201323038
  * Andrews et al. (2017) Andrews, B. H., Weinberg, D. H., Schönrich, R., & Johnson, J. A. 2017, ApJ, 835, 224, doi: 10.3847/1538-4357/835/2/224
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Barbuy et al. (2018) Barbuy, B., Chiappini, C., & Gerhard, O. 2018, ARA&A, 56, 223, doi: 10.1146/annurev-astro-081817-051826
  * Barnes (2007) Barnes, S. A. 2007, ApJ, 669, 1167, doi: 10.1086/519295
  * Bashi & Zucker (2019) Bashi, D., & Zucker, S. 2019, AJ, 158, 61, doi: 10.3847/1538-3881/ab27c9
  * Bashi et al. (2020) Bashi, D., Zucker, S., Adibekyan, V., et al. 2020, A&A, 643, A106, doi: 10.1051/0004-6361/202038881
  * Belokurov et al. (2020) Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880, doi: 10.1093/mnras/staa876
  * Bennett & Bovy (2019) Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417, doi: 10.1093/mnras/sty2813
  * Bensby & Feltzing (2006) Bensby, T., & Feltzing, S. 2006, MNRAS, 367, 1181, doi: 10.1111/j.1365-2966.2006.10037.x
  * Bensby et al. (2003) Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527, doi: 10.1051/0004-6361:20031213
  * Bensby et al. (2014) Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71, doi: 10.1051/0004-6361/201322631
  * Bensby et al. (2007) Bensby, T., Zenn, A. R., Oey, M. S., & Feltzing, S. 2007, ApJ, 663, L13, doi: 10.1086/519792
  * Bird et al. (2013) Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43, doi: 10.1088/0004-637X/773/1/43
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/aa7567
  * Bond et al. (2010) Bond, N. A., Ivezić, Ž., Sesar, B., et al. 2010, ApJ, 716, 1, doi: 10.1088/0004-637X/716/1/1
  * Bournaud et al. (2009) Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJ, 707, L1, doi: 10.1088/0004-637X/707/1/L1
  * Bovy et al. (2012) Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148, doi: 10.1088/0004-637X/753/2/148
  * Brook et al. (2006) Brook, C. B., Kawata, D., Martel, H., Gibson, B. K., & Bailin, J. 2006, ApJ, 639, 126, doi: 10.1086/499154
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765, doi: 10.1086/303726
  * Ciucă et al. (2021) Ciucă, I., Kawata, D., Miglio, A., Davies, G. R., & Grand, R. J. J. 2021, MNRAS, 503, 2814, doi: 10.1093/mnras/stab639
  * Connolly et al. (2000) Connolly, A. J., Genovese, C., Moore, A. W., et al. 2000, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0008187
  * Duong et al. (2018) Duong, L., Freeman, K. C., Asplund, M., et al. 2018, MNRAS, 476, 5216, doi: 10.1093/mnras/sty525
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Freeman & Bland-Hawthorn (2002) Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487, doi: 10.1146/annurev.astro.40.060401.093840
  * Fuhrmann (2004) Fuhrmann, K. 2004, Astronomische Nachrichten, 325, 3, doi: 10.1002/asna.200310173
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gilmore & Reid (1983) Gilmore, G., & Reid, N. 1983, MNRAS, 1025, doi: 10.1093/mnras/202.4.1025
  * Gravity Collaboration et al. (2018) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15, doi: 10.1051/0004-6361/201833718
  * Grisoni et al. (2017) Grisoni, V., Spitoni, E., Matteucci, F., et al. 2017, MNRAS, 472, 3637, doi: 10.1093/mnras/stx2201
  * 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
  * Haywood et al. (2013) Haywood, M., Di Matteo, P., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109, doi: 10.1051/0004-6361/201321397
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Isobe et al. (1986) Isobe, T., Feigelson, E. D., & Nelson, P. I. 1986, ApJ, 306, 490, doi: 10.1086/164359
  * Jørgensen & Lindegren (2005) Jørgensen, B. R., & Lindegren, L. 2005, A&A, 436, 127, doi: 10.1051/0004-6361:20042185
  * Katz et al. (2011) Katz, D., Soubiran, C., Cayrel, R., et al. 2011, A&A, 525, A90, doi: 10.1051/0004-6361/201014840
  * Kawata et al. (2018) Kawata, D., Allende Prieto, C., Brook, C. B., et al. 2018, MNRAS, 473, 867, doi: 10.1093/mnras/stx2464
  * Kazantzidis et al. (2008) Kazantzidis, S., Bullock, J. S., Zentner, A. R., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254, doi: 10.1086/591958
  * Kordopatis et al. (2011) Kordopatis, G., Recio-Blanco, A., de Laverny, P., et al. 2011, A&A, 535, A107, doi: 10.1051/0004-6361/201117373
  * Lagarde et al. (2021) Lagarde, N., Reylé, C., Chiappini, C., et al. 2021, A&A, 654, A13, doi: 10.1051/0004-6361/202039982
  * Lee et al. (2011) Lee, Y. S., Beers, T. C., An, D., et al. 2011, ApJ, 738, 187, doi: 10.1088/0004-637X/738/2/187
  * Leung & Bovy (2019) Leung, H. W., & Bovy, J. 2019, MNRAS, 489, 2079, doi: 10.1093/mnras/stz2245
  * Lian et al. (2020) Lian, J., Thomas, D., Maraston, C., et al. 2020, MNRAS, 497, 2371, doi: 10.1093/mnras/staa2078
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Mackereth et al. (2019) Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176, doi: 10.1093/mnras/stz1521
  * Marín-Franch et al. (2009) Marín-Franch, A., Aparicio, A., Piotto, G., et al. 2009, ApJ, 694, 1498, doi: 10.1088/0004-637X/694/2/1498
  * Miglio et al. (2021) Miglio, A., Chiappini, C., Mackereth, J. T., et al. 2021, A&A, 645, A85, doi: 10.1051/0004-6361/202038307
  * Minchev et al. (2015) Minchev, I., Martig, M., Streich, D., et al. 2015, ApJ, 804, L9, doi: 10.1088/2041-8205/804/1/L9
  * Minchev et al. (2019) Minchev, I., Matijevic, G., Hogg, D. W., et al. 2019, MNRAS, 487, 3946, doi: 10.1093/mnras/stz1239
  * Ness et al. (2013) Ness, M., Freeman, K., Athanassoula, E., et al. 2013, MNRAS, 430, 836, doi: 10.1093/mnras/sts629
  * Nikakhtar et al. (2021) Nikakhtar, F., Sanderson, R. E., Wetzel, A., et al. 2021, ApJ, 921, 106, doi: 10.3847/1538-4357/ac1a10
  * Nissen & Schuster (2010) Nissen, P. E., & Schuster, W. J. 2010, A&A, 511, L10, doi: 10.1051/0004-6361/200913877
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science and Engineering, 9, 10, doi: 10.1109/MCSE.2007.58
  * Pinsonneault et al. (2018) Pinsonneault, M. H., Elsworth, Y. P., Tayar, J., et al. 2018, ApJS, 239, 32, doi: 10.3847/1538-4365/aaebfd
  * Re Fiorentin et al. (2019) Re Fiorentin, P., Lattanzi, M. G., & Spagna, A. 2019, MNRAS, 484, L69, doi: 10.1093/mnrasl/sly198
  * 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
  * Reddy et al. (2006) Reddy, B. E., Lambert, D. L., & Allende Prieto, C. 2006, MNRAS, 367, 1329, doi: 10.1111/j.1365-2966.2006.10148.x
  * Rojas-Arriagada et al. (2016) Rojas-Arriagada, A., Recio-Blanco, A., de Laverny, P., et al. 2016, A&A, 586, A39, doi: 10.1051/0004-6361/201526969
  * Ruchti et al. (2010) Ruchti, G. R., Fulbright, J. P., Wyse, R. F. G., et al. 2010, ApJ, 721, L92, doi: 10.1088/2041-8205/721/2/L92
  * Sanders & Das (2018) Sanders, J. L., & Das, P. 2018, MNRAS, 481, 4093, doi: 10.1093/mnras/sty2490
  * Schönrich & Binney (2009a) Schönrich, R., & Binney, J. 2009a, MNRAS, 399, 1145, doi: 10.1111/j.1365-2966.2009.15365.x
  * Schönrich & Binney (2009b) —. 2009b, MNRAS, 396, 203, doi: 10.1111/j.1365-2966.2009.14750.x
  * 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
  * Schuster et al. (2012) Schuster, W. J., Moreno, E., Nissen, P. E., & Pichardo, B. 2012, A&A, 538, A21, doi: 10.1051/0004-6361/201118035
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461
  * Shao & Zhao (1996) Shao, Z., & Zhao, J. 1996, Acta Astronomica Sinica, 37, 377
  * Sharma et al. (2021) Sharma, S., Hayden, M. R., & Bland-Hawthorn, J. 2021, MNRAS, 507, 5882, doi: 10.1093/mnras/stab2015
  * Soderblom (2010) Soderblom, D. R. 2010, ARA&A, 48, 581, doi: 10.1146/annurev-astro-081309-130806
  * Soltis et al. (2021) Soltis, J., Casertano, S., & Riess, A. G. 2021, ApJ, 908, L5, doi: 10.3847/2041-8213/abdbad
  * Soubiran & Girard (2005) Soubiran, C., & Girard, P. 2005, A&A, 438, 139, doi: 10.1051/0004-6361:20042390
  * Spitoni et al. (2019) Spitoni, E., Silva Aguirre, V., Matteucci, F., Calura, F., & Grisoni, V. 2019, A&A, 623, A60, doi: 10.1051/0004-6361/201834188
  * Spitoni et al. (2020) Spitoni, E., Verma, K., Silva Aguirre, V., & Calura, F. 2020, A&A, 635, A58, doi: 10.1051/0004-6361/201937275
  * Spitoni et al. (2021) Spitoni, E., Verma, K., Silva Aguirre, V., et al. 2021, A&A, 647, A73, doi: 10.1051/0004-6361/202039864
  * Takeda et al. (2007) Takeda, G., Ford, E. B., Sills, A., et al. 2007, ApJS, 168, 297, doi: 10.1086/509763
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vincenzo et al. (2021) Vincenzo, F., Weinberg, D. H., Miglio, A., Lane, R. R., & Roman-Lopes, A. 2021, MNRAS, 508, 5903, doi: 10.1093/mnras/stab2899
  * Weinberg et al. (2017) Weinberg, D. H., Andrews, B. H., & Freudenburg, J. 2017, ApJ, 837, 183, doi: 10.3847/1538-4357/837/2/183
  * Wojno et al. (2016) Wojno, J., Kordopatis, G., Steinmetz, M., et al. 2016, MNRAS, 461, 4246, doi: 10.1093/mnras/stw1633
  * Wyse & Gilmore (1988) Wyse, R. F. G., & Gilmore, G. 1988, AJ, 95, 1404, doi: 10.1086/114738
  * Yu et al. (2021) Yu, Z., Li, J., Chen, B., et al. 2021, ApJ, 912, 106, doi: 10.3847/1538-4357/abf098