# The Local Dark Matter Kinematic Substructure Based on LAMOST K Giants

Hai Zhu  School of Physics and Astronomy, Shanghai Jiao Tong University, Key
Laboratory for Particle Astrophysics and Cosmology (MoE), Shanghai Key
Laboratory for Particle Physics and Cosmology, Shanghai 200240, China  Rui Guo
School of Physics and Astronomy, Shanghai Jiao Tong University, Key Laboratory
for Particle Astrophysics and Cosmology (MoE), Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai 200240, China  Juntai Shen  School of
Physics and Astronomy, Shanghai Jiao Tong University, Key Laboratory for
Particle Astrophysics and Cosmology (MoE), Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai 200240, China  Jianglai Liu  New
Cornerstone Science Laboratory, Tsung-Dao Lee Institute, Shanghai Jiao Tong
University, Shanghai, 200240, China  School of Physics and Astronomy, Shanghai
Jiao Tong University, Key Laboratory for Particle Astrophysics and Cosmology
(MoE), Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai
200240, China  Shanghai Jiao Tong University Sichuan Research Institute,
Chengdu 610213, China  Chao Liu  Key Laboratory of Space Astronomy and
Technology, National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100101, China  Institute for Frontiers in Astronomy and Astrophysics,
Beijing Normal University, Beijing 102206, China  Xiang-Xiang Xue  CAS Key
Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese
Academy of Sciences, Beijing 100101, China  Institute for Frontiers in
Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, China
Lan Zhang  CAS Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China  Shude Mao
Department of Astronomy and Tsinghua Centre for Astrophysics, Tsinghua
University, Beijing 100084, China

###### Abstract

Numerical simulations indicate that correlations exist between the velocity
distributions of stars and dark matter (DM). We study the local DM velocity
distribution based on these correlations. We select K giants from LAMOST DR8
cross-matched with Gaia DR3, which has robust measurements of three-
dimensional velocity and metallicity, and separate them into the disk, halo
substructure and main halo components in the chemo-dynamical space utilizing
the Gaussian Mixture Model. The substructure component is highly radially
anisotropic, and possibly related to the Gaia-Enceladus-Sausage (GES) merger
event, while the halo component is isotropic and accreted from the earliest
mergers following the Maxwell-Boltzmann Distribution (Standard Halo Model,
SHM). We find that the GES-like substructure contributes $\sim 85\%$ of the
local non-disk stars in the Solar neighbourhood, which is nearly invariant
when applying different volume cuts or additional angular momentum
constraints. Utilizing the metallicity-stellar-mass relation and the stellar-
mass-halo-mass relation, we find that $\sim 25_{-15}^{+24}\%$ of local DM is
in the kinematic substructure. Combined with the stellar distributions of non-
disk components, we compute the velocity distribution of local DM. The
modified heliocentric velocity distribution of local DM shifts to a lower
speed and has a sharper peak compared to the SHM, which yields updated
detection limits for the DM direct detection experiments. Our work confirms
that the local DM velocity distribution deviates from the SHM, and needs to be
properly accounted in the DM detection experiments.

Galaxy: stellar halos – Galaxy: kinematics and dynamics

##  1 Introduction

In the $\Lambda$CDM paradigm, the Milky Way halo is built up from hierarchical
merger of smaller satellites (White & Rees, 1978). Depending on the merger
time and properties of the satellites, such as the mass ratio and the
satellite orbits, the satellite remnants that construct the Milky Way halo at
present day can be separated into different types. The earliest merger has the
longest time to virialize after the dynamical relaxation, which results in a
smooth and isotropic distribution. This isotropic component can be described
by the Standard Halo Model (SHM) (Ostriker et al., 1974; Bahcall & Soneira,
1980; Caldwell & Ostriker, 1981) that follows the Maxwell-Boltzmann velocity
distribution (Drukier et al., 1986; Freese et al., 1988). The more recent
merger events usually do not have enough time to fully mix into the main halo
even in the configuration space. In particular, the outer part of the
satellite is often stripped away to form a tidal stream. The stellar streams
have been widely used to study the Milky Way mass distribution and its dark
matter (DM) halo shape (e.g. Malhan & Ibata, 2019; Vasiliev et al., 2020;
Hattori et al., 2021; Ibata et al., 2021). An intermediate merger between
ancient and recent leaves substructures called the “debris flows”, which are
well mixed with the main halo in the configuration space, but showing
structures in the kinematic space (Helmi & White, 1999; Gómez et al., 2010).

The most important merger event that built up the Milky Way might be the Gaia-
Enceladus-Sausage (GES) (Helmi et al., 2018; Belokurov et al., 2018), which
remains a sausage-like structure in the radial and azimuthal
$v_{r}$-$v_{\phi}$ velocity space. The progenitor of the GES is thought to be
a large dwarf galaxy with a highly eccentric orbit and a stellar mass on the
order of $10^{9}-10^{10}$ M⊙ (Vincenzo et al., 2019; Deason et al., 2019;
Fattahi et al., 2019; Koppelman et al., 2020) and merged with the Milky Way
around 8$-$10 Gyr ago (Sahlholdt et al., 2019; Bignone et al., 2019; Bonaca et
al., 2020). The chemo-dynamical studies indicate that the GES stars have a
metallicity peak of [Fe/H] around $-1.4$ to $-1.2$ dex (Sahlholdt et al.,
2019; Gallart et al., 2019; Das et al., 2020; Feuillet et al., 2020; Naidu et
al., 2020). Due to its large mass, GES strongly influenced the radial
kinematic profiles of the stellar halo. This highly radially anisotropic
component exists in the local stellar halo (within $\sim$ 10 kpc from the Sun)
and extends out to a Galactocentric radius $r_{gc}\sim$ 25 kpc (Lancaster et
al., 2019; Bird et al., 2019, 2021; Wu et al., 2022). Deason et al. (2018)
showed that the apocenter of the GES debris stars is about $20-25~{}$kpc,
which might be related to the break radius of the density profile of the
Galactic stellar halo (Watkins et al., 2009; Sesar et al., 2011; Pila-Díez, B.
et al., 2015; Xue et al., 2015). Wu et al. (2022) found that the GES
contributes about $41\%-74\%$ of the inner ($r_{gc}<30$ kpc) stellar halo. In
the Solar neighbourhood, Necib et al. (2019a) showed that the GES stars
contribute $60\%-80\%$ of the non-disk components of main sequence stars in
the SDSS DR9 sample, with a peak metallicity of $\rm{[Fe/H]}=-1.4$ dex.

In the DM direct detection experiments, both the density and velocity
distributions of the local DM are important inputs to calculate the
differential scattering rates and the detection limit curves (e.g. Lewin &
Smith, 1996; Frandsen et al., 2012; Bhattacharjee et al., 2013; Fairbairn et
al., 2013; Green, 2017). The local DM density has been widely measured
utilizing the local stellar kinematics and the Galactic rotation curve (e.g.
Xia et al., 2016; Hagen, Jorrit H. J. & Helmi, Amina, 2018; Huang et al.,
2016; de Salas et al., 2019; Eilers et al., 2019; Guo et al., 2020, 2022),
while the DM velocity distribution is usually assumed to follow the Maxwell-
Boltzmann distribution. However, the prior merger events and the recent
interactions with satellites could change the local DM velocity distribution
(Necib et al., 2019a; Donlon et al., 2022; Besla et al., 2019; Smith-Orlik et
al., 2023; Donaldson et al., 2022; Evans et al., 2019). Evans et al. (2019)
found that the strongly radially anisotropic GES may increase the fraction of
high velocity tail relative to the SHM, leading to a harder kinetic energy
distribution of incoming DM particles in detection experiments and hence a
strong exclusion limit for lower mass DM. Nevertheless, Necib et al. (2019a)
utilized the stellar velocity distribution to trace the DM velocity
distribution, and found a contrary influence of the GES on the DM detection
experiments. The differences mainly result from the contribution of the GES to
the local DM and the profiles of its velocity distribution. This tension
requests for a more detailed study on the exact stellar velocity distribution
of the halo and substructures.

Some simulations of Milky-Way-like galaxies show that there are correlations
between the velocity distribution of stars and the accreted DM particles
(Herzog-Arbeitman et al., 2018a; Lisanti et al., 2015; Necib et al., 2019b),
i.e. the stars and DM accreted from the same merger satellite may share
similar velocity distribution111The stars formed in the accreted gas are not
included, as their kinematic information could be different after the star
formation.. This allows one to trace the DM velocity distribution with halo
stars. The velocity distribution of the DM accreted from the oldest mergers
forming the virialized halo can be well traced by metal-poor stars in the Eris
hydrodynamic simulations (Herzog-Arbeitman et al., 2018a), which motivates the
first study using the RAVE-TGAS data to recover the velocity distribution of
the local relaxed DM component in Herzog-Arbeitman et al. (2018b). The DM
accreted from younger mergers, which remained as debris flows and streams, are
associated with more metal-rich stars, as illustrated in the Via Lactea
simulations and FIRE-2 simulations (Lisanti et al., 2015; Necib et al.,
2019b). The correlation between DM and star particles found in the FIRE-2
simulations holds for two simulated hosts with different merger histories,
motivating the study of the kinematic substructure of DM in the Solar
neighbourhood in Necib et al. (2019a). In this framework the stellar
metallicity distribution is crucial and requires a large sample of stars with
accurate metallicity measurement. Here we attempt to provide an independent
check on the contribution of the GES to the stellar halo utilizing another
stellar tracer, K giants selected from LAMOST (Cui et al., 2012; Zhao et al.,
2012). These K giants have high-quality measurements of line-of-sight velocity
and metallicity, and are cross-matched with Gaia DR3 (Gaia Collaboration et
al., 2022) for proper motions measurements. More importantly, our K giants
sample has accurate distance estimation. We follow the study of Necib et al.
(2019a) to extract information of different components in the chemo-dynamical
space applying the Gaussian Mixture Model (GMM) method. Assuming that there
exists correlations between the DM and halo stars, we obtain the velocity
distribution of the local DM considering the kinematic substructure of the
GES, and then to investigate its influence on the DM direct detection
experiments.

This paper is organized as follows: In Section 2, we discuss the selection of
the K giants sample and the likelihood calculation procedure of the GMM
method. In Section 3, we present the distributions of different stellar
components from our best-fit model and explore the variation due to different
spatial volume cuts of the samples. To reduce the dominance of the disk stars,
we apply an angular momentum cut to the samples, and extend the GMM fitting to
halo stars near the disk in Section 4. The application of the results to the
DM velocity distribution and direct detection experiments are shown in Section
5. Discussions and summary are presented in Section 6 and Section 7,
respectively.

##  2 Data Analysis

###  2.1 Sample Selection

In order to study the kinematic substructure in the local stellar halo, we
consider the K giants selected from LAMOST Data Release 8 (DR8). The selection
criteria of the LAMOST K giants are based on the surface gravity log $g$ and
the effective temperature $T_{\rm{eff}}$ according to Liu et al. (2014), i.e.
log $g$ $\leq 4$ when $4600\leq$$T_{\rm{eff}}$/K$\leq 5600$ or log $g$ $\leq
3.5$ when 4000 $\leq$$T_{\rm{eff}}$/K$\leq 4600$. We cross-match the K giant
sample with Gaia DR3 with a radius $1^{\prime\prime}$ in order to obtain
accurate kinematic information. LAMOST provides accurate metallicity ([Fe/H])
and line-of-sight velocity ($v_{\rm{los}}$) measurements from the spectra,
while Gaia provides accurate proper motions. The cross-matched sample contains
612,136 K giants. The LAMOST $v_{\rm{los}}$ has a systematic offset of 5.34 km
s-1 for K giants relative to that from Gaia, which was investigated (Tian et
al., 2015; Ding et al., 2021) and has been corrected in this work. Due to the
large uncertainty in distance derived from parallax for distant stars, the
distance ($d$) is estimated with the Bayesian method described in Xue et al.
(2014). We further exclude stars with distance error $\delta d/d>0.3$ and
absolute spherical velocity $|\boldsymbol{v}|>400$ km s-1. This results in a
sample containing 438,466 K giants, 18,484 of which are potential halo stars
within the Solar neighbourhood ($r\in[7.5,10]$ kpc & $|z|>2.5$ kpc). This K
giant sample has a mean uncertaintiy of 20% in the distance, 4.4 km s-1 in
$v_{\rm{los}}$ and 0.06 dex in the metallicity. The metallicity is cut down to
$-2.5$ dex due to the limit of LAMOST spectroscopic pipelines (Wu et al.,
2011; Luo et al., 2015).

The coordinates system adopted in this paper is the Galactocentric spherical
coordinates ($r,\theta,\phi$), pointing outwards, to the north Galactic pole
and to the Galactic rotation, respectively. We adopt a Solar distance to the
Galactic center of $R_{0}=8.122$ kpc (GRAVITY Collaboration et al., 2018), and
a Solar vertical height of $z_{\odot}=20.8$ pc (Bennett & Bovy, 2018), with a
Solar motion of ($+12.9,+245.6,+7.78$) km s-1 in the Galactocentric Cartesian
coordinates (Reid & Brunthaler, 2004; GRAVITY Collaboration et al., 2018;
Drimmel & Poggio, 2018). We utilize the python package astropy (Astropy
Collaboration et al., 2018) to transform the astrometric measurements of
proper motions ($\mu_{\alpha},\mu_{\delta}$) and $v_{\rm{los}}$ to the
kinematic information in the spherical coordinates. The measurement errors are
propagated to the uncertainties in the velocities using a Monte Carlo sampling
with the python package pyia (Price-Whelan, 2018). Fig. 1 shows the spatial
distribution of the sample we used. We focus mainly on the region of
$r\in[7.5,10.0]$ kpc and $|z|>$2.5 kpc, and explore the regions closer to the
Galactic plane in Section 4.

Figure 1: Spatial distribution of the LAMOST DR8 $\times$ Gaia DR3 sample in
terms of Galactocentric radius $r$ and vertical distance from the Galactic
disk $z$. Stars within $r\in[7.5,10.0]$ kpc and $|z|>$ 2.5 kpc are utilized in
the following chemo-dynamical analyses.

###  2.2 Procedure of Likelihood Calculation

The stellar sample can be separated into three components, i.e. the halo, the
substructure and the disk, according to the stellar chemo-dynamical
information (Necib et al., 2019a). The halo component is referred to the
isotropic and metal-poor halo contributed by the oldest mergers, while the
substructure is a radially anisotropic halo component mainly contributed by
the GES merger event. The $v_{r}$ distribution of the substructure shows a
symmetric bimodal Gaussian distribution, due to the highly radial merger
trajectory. The disk component contains more metal-rich stars born in-situ
with strong azimuthal rotation. As these components show different chemical
and kinematic properties, we utilize the GMM to identify different components
with a Bayesian estimation (Lancaster et al., 2019; Necib et al., 2019a), in
the ($v_{r},\ v_{\theta},\ v_{\phi},\ \rm\rm{[Fe/H]}$) chemo-dynamical space.

We follow the framework in Necib et al. (2019a) to define the likelihood of
each star $\mathcal{L}_{i}$, associated with a set of observables labeled as
$O_{i}=(\boldsymbol{v}_{i},~{}{\rm{[Fe/H]}}_{i})$ and different model
parameters. We label the disk, halo and substructure components with flag of
$j=d,h,s$, respectively.

The likelihood of a star belonging to the disk component is

$\displaystyle\mathcal{L}_{d}(O_{i}|\Theta)$ (1)
$\displaystyle=\mathcal{N}(\boldsymbol{v}_{i}|\boldsymbol{\mu}^{d},\boldsymbol{\Sigma}_{i}^{d})\mathcal{N}({\rm{[Fe/H]}}_{i}|\mu^{d}_{{\rm{[Fe/H]}}},\sigma_{{\rm{[Fe/H]}},i}^{d})\
.$

Similarly, the likelihood function of the halo component is

$\displaystyle\mathcal{L}_{h}(O_{i}|\Theta)$ (2)
$\displaystyle=\mathcal{N}(\boldsymbol{v}_{i}|\boldsymbol{\mu}^{h},\boldsymbol{\Sigma}_{i}^{h})\mathcal{N}({\rm{[Fe/H]}}_{i}|\mu^{h}_{{\rm{[Fe/H]}}},\sigma_{{\rm{[Fe/H]}},i}^{h})\
,$

where $\Theta$ represents the parameter set and $\mathcal{N}$ denotes the
Gaussian distribution (3 dimensional Gaussian distribution for the velocity).
$\Theta$ contains the free parameters of the Gaussian distributions for both
the velocity and metallicity. ${\boldsymbol{\mu}^{j}}$ is the mean velocity
($\mu_{r,j},\mu_{\theta,j},\mu_{\phi,j}$), while $\boldsymbol{\Sigma}^{j}$ is
the covariance matrix determined by the velocity dispersion
$\sigma_{r,\theta,\phi}$ and the correlation coefficients
$\rho_{r\theta,r\phi,\theta\phi}$. The correlation coefficients are found to
be very small (e.g. Evans et al., 2016; Necib et al., 2019a; Wu et al., 2022).
In observations, the dispersion of metallicity and the velocity covariance
matrix is a combination of the true value and the measurement errors, i.e.
$\boldsymbol{\Sigma}_{\rm{obs}}=\boldsymbol{\Sigma}_{\rm{true}}+\boldsymbol{\Sigma}_{\rm{err}},\sigma_{{\rm{[Fe/H]}},\rm{obs}}=\sqrt{\sigma_{{\rm{[Fe/H]}},\rm{true}}^{2}+\sigma_{{\rm{[Fe/H]}},\rm{err}}^{2}}$,
varying for each star. Considering the mean and covariance matrix of velocity,
and the mean and dispersion of metallicity, the models of disk and halo
components each contain 11 free parameters.

Belokurov et al. (2018) and Necib et al. (2019a) demonstrated that the
velocity of the substructure component resulted from the GES is highly
radially anisotropic. Thus it can be described by a bimodal Gaussian
distribution with a likelihood of

$\displaystyle\mathcal{L}_{s}(O_{i}|\Theta)=$
$\displaystyle\frac{1}{2}[\mathcal{N}(\boldsymbol{v}_{i}|\boldsymbol{\mu}^{\tilde{s}},\boldsymbol{\Sigma}_{i}^{s})+\mathcal{N}(\boldsymbol{v}_{i}|\boldsymbol{\mu}^{s},\boldsymbol{\Sigma}_{i}^{s})]$
(3)
$\displaystyle\times\mathcal{N}({\rm{[Fe/H]}}_{i}|\mu^{s}_{{\rm{[Fe/H]}}},\sigma_{{\rm{[Fe/H]}},i}^{s})\
,$

where $\boldsymbol{\mu}^{s}\\!=\\!(\mu_{r,s},\mu_{\theta,s},\mu_{\phi,s})$ and
$\boldsymbol{\mu}^{\tilde{s}}\\!=\\!(-\mu_{r,s},\mu_{\theta,s},\mu_{\phi,s})$.
This model resembles a single Gaussian distribution in $v_{r}$ when
$\mu_{r,s}$ is small and $\sigma_{r,s}$ is large, in which case it results in
a strong degeneracy between the halo and the substructure components. The
model of the substructure component gives another 11 free parameters.

Table 1: Priors of 35 free parameters. Parameter | Type | Priors
---|---|---
Disk | Halo | Substructure
$\mu_{r}$ | Linear | $[-70,\ 70]$ | $[-70,\ 70]$ | $[0,\ 250]$
$\mu_{\theta}$ | Linear | $[-70,\ 70]$ | $[-70,\ 70]$ | $[-70,\ 70]$
$\mu_{\phi}$ | Linear | $[0,\ 300]$ | $[-70,\ 70]$ | $[-70,\ 70]$
$\sigma_{r,\theta,\phi}$ | Linear | $[0,\ 200]$ | $[0,\ 200]$ | $[0,\ 200]$
$\rho_{r\theta,r\phi,\theta\phi}$ | Linear | $[-1,\ 1]$ | $[-1,\ 1]$ | $[-1,\ 1]$
$\mu_{\mathrm{[Fe/H]}}$ | Linear | $[-1,\ 0]$ | $[-2.5,\ -1.5]$ | $[-1.5,\ -1]$
$\sigma_{\mathrm{[Fe/H]}}$ | Linear | $[0,\ 2]$ | $[0,\ 2]$ | $[0,\ 2]$
$\mathcal{Q}$ | Linear | $[0,\ 1]$ | $\cdots$ | $[0,\ 1]$

Thus, the likelihood for a sample containing $N$ stars is defined as

$\mathcal{L}(\\{O_{i}\\}|\Theta)=\prod_{i=1}^{N}\sum_{j=d,h,s}\mathcal{Q}_{j}\mathcal{L}_{j}(O_{i}|\Theta)\
,$ (4)

where $\mathcal{Q}_{j}$ denotes the probability that the star belongs to the
j-th component with an normalization constraint of
$\mathcal{Q}_{h}=1-\mathcal{Q}_{d}-\mathcal{Q}_{s}$. $\mathcal{Q}_{j}$ also
represents the integrated fractional area of any marginalized 1D distribution
for that component. With the two fraction parameters of $\mathcal{Q}_{d}$ and
$\mathcal{Q}_{s}$, we have 35 free parameters in total.

We apply the emcee (Foreman-Mackey et al., 2013) package to perform the Markov
chain Monte Carlo (MCMC) Bayesian estimation to find the posterior
distributions of all 35 free parameters. We utilize 250 walkers and 5000 steps
as burn-in, followed by another 5000 steps to obtain the posterior
distributions of parameters. The flat priors of each parameter are listed in
Table 1, similar to Necib et al. (2019a). As the halo component contributes
less in the sample, we utilize narrower prior ranges for the mean
metallicities based on previous results (Necib et al., 2019a; Wu et al.,
2022). We adopt the set of parameters with the maximum likelihood as the best-
fit result and the 2.1%, 50%, 97.9% percentiles of the 1D marginalized
posterior distributions to denote the median and $\sim 2\sigma$ uncertainties.

Figure 2: Radial velocity ($v_{r}$) vs. azimuthal velocity ($v_{\phi}$) color-
coded by [Fe/H] for the stellar samples within $r\in[7.5,10]$ kpc $\&\
|z|>2.5$ kpc (left) and $r\in[7.5,10]$ kpc $\&\ |z|>5.0$ kpc (right).

##  3 The Stellar Distribution in the Chemo-dynamical Space

As illustrated by Schönrich & Binney (2011) and Necib et al. (2019a), below
$|z|\sim 2.5$ kpc, the azimuthal velocity distribution of the disk cannot be
fitted with a single Gaussian distribution due to the contribution from the
thin disk. Thus, in this Section, we focus on the analyses of the stars within
$r\in[7.5,10.0]$ kpc and $|z|>$ 2.5 kpc. The velocity distribution of $v_{r}$
versus $v_{\phi}$ of this sample is displayed in the left panel of Fig. 2. The
clump with a high metallicity and an azimuthal rotation of $\sim 200$ km s-1
is the disk component, which is the dominant component in this sample. The
elongated “sausage” structure in $v_{r}$ direction with a lower metallicity is
the substructure component contributed by the GES. It becomes the dominant
component when a larger vertical cut of $|z|>$ 5 kpc is applied (right panel
of Fig. 2) where the disk component disappears. In order to extrapolate the
result to the Galactic plane, we apply an additional angular momentum
($L_{z}$) cut to reduce the contribution of the disk component in the later
Section 4.

###  3.1 The Solar Neighbourhood

Figure 3: The best-fit chemo-dynamical distribution
($v_{r},v_{\theta},v_{\phi},\rm{[Fe/H]}$) of the sample in the Solar
neighbourhood with $r\in[7.5,8.5]~{}$kpc & $|z|>2.5~{}$kpc. The gray histogram
is the data, while the colored lines are from the model which represent the
disk, substructure and halo components as green, blue and red lines,
respectively.

As we aim to obtain the local DM velocity distribution, we select a sample in
the Solar neighbourhood with a radial range of $r\in[7.5,8.5]~{}$kpc and a
vertical cut of $|z|>2.5~{}$kpc. We apply the GMM method and the MCMC
technique to separate the three components. The results are shown in Fig. 3,
and the posterior distributions of the free parameters are displayed in the
Appendix for each component separately. Our best-fit model well fits the
chemo-dynamical distributions of the sample. The disk component is the most
dominant, while the halo component contributes the least. The fractions of the
three components are $\mathcal{Q}_{d}=0.55^{+0.02}_{-0.02}$,
$\mathcal{Q}_{s}=0.39^{+0.02}_{-0.02}$ and
$\mathcal{Q}_{h}=0.06^{+0.02}_{-0.01}$, for the disk, substructure and halo
components, respectively. The relative fraction between the substructure and
halo components is
$\mathcal{Q}_{s}/(\mathcal{Q}_{h}+\mathcal{Q}_{s})=0.87^{+0.03}_{-0.04}$,
which means that the substructure contributes $\sim 87\%$ in all local halo
stars. It should be noted that there are potential impacts of the selection
effect of the spectroscopic survey (sample incompleteness compared to the
photometric survey) on this fraction. Wu et al. (2022) shows that the LAMOST
selection function mainly affects the GES fraction
$\mathcal{Q}_{s}/(\mathcal{Q}_{h}+\mathcal{Q}_{s})$ for fainter stars in the
very outer halo ($r<30$ kpc). For brighter stars in the Solar neighbourhood,
the impact of selection function is negligible as shown in their Fig. 24. The
selection effect are not taken into consideration in this study since our K
giants are within $r<10$ kpc.

The disk component is relatively metal-rich with a mean metallicity of
$\mu_{\rm{[Fe/H]}}=-0.55_{-0.01}^{+0.01}$ dex and a dispersion of
$\sigma_{\rm{[Fe/H]}}=0.17_{-0.01}^{+0.01}$ dex. Its motion is dominated by
the azimuthal rotation with $\mu_{\phi}=134.73_{-2.87}^{+2.74}$ km s-1 and
$\sigma_{\phi}=65.89_{-2.17}^{+2.12}$ km s-1. The substructure component is
radially anisotropic, showing a bimodal Gaussian distribution in $v_{r}$ with
$\mu_{r}=120.09_{-6.61}^{+6.44}$ km s-1 and
$\sigma_{r}=108.52_{-4.53}^{+5.20}$ km s-1. As we adopt the right-handed
Galactocentric frame, $v_{\phi}>0$ means a prograde rotation. Thus, for the
substructure component, $\mu_{\phi}=22.59_{-4.61}^{+4.74}$ km s-1 and
$\sigma_{\phi}=93.43_{-3.29}^{+3.20}$ km s-1 indicates a slightly prograde
orbit as referred in some previous works about the debris flow of the massive
radial merger event GES (Belokurov et al., 2018; Necib et al., 2019a; Wu et
al., 2022), whereas other works show a retrograde rotation (e.g. Helmi et al.,
2018; Koppelman et al., 2020; Mackereth et al., 2018). Our sample predicts a
more metal-rich substructure component with
$\mu_{\rm{[Fe/H]}}=-1.15_{-0.04}^{+0.04}$ dex than the previous result of
$\mu_{\rm{[Fe/H]}}=-1.39$ dex (Necib et al., 2019a). The difference could
result from the degeneracy between the substructure and halo components. The
halo component consists of the accreted stars from the oldest mergers, which
thus becomes relaxed and isotropic with almost zero mean velocities and large
dispersions, i.e.
$(\mu_{r},\mu_{\theta},\mu_{\phi})=(-4.59_{-18.19}^{+16.53},8.02_{-13.30}^{+13.74},19.82_{-14.06}^{+14.19})$
km s-1 and $(\sigma_{r},\sigma_{\theta},\\\
\sigma_{\phi})=(128.47_{-12.97}^{+12.85},113.62_{-9.52}^{+11.59},112.02_{-9.69}^{+11.00})$
km s-1. Note that the halo component also shows a slightly prograde azimuthal
rotation similar to the substructure component. The metallicity of the halo
component is more metal-poor than many previous works with
$\mu_{\rm{[Fe/H]}}=-2.04_{-0.06}^{+0.07}$ dex and
$\sigma_{\rm{[Fe/H]}}=0.21_{-0.04}^{+0.04}$ dex.

All of the velocity correlation coefficients for each component are small as
expected: for the disk component
$(\rho_{r\theta},\rho_{r\phi},\rho_{\theta\phi})=(0.16_{-0.04}^{+0.04},-0.09_{-0.04}^{+0.04},0.02_{-0.04}^{+0.04})$;
for the substructure component
$(\rho_{r\theta},\rho_{r\phi},\rho_{\theta\phi})\\!\\!=\\!\\!(0.15_{-0.05}^{+0.05},-0.08_{-0.06}^{+0.06},0.06_{-0.04}^{+0.04})$;
for the halo component
$(\rho_{r\theta},\rho_{r\phi},\rho_{\theta\phi})=(-0.07_{-0.12}^{+0.12},0.03_{-0.13}^{+0.16},-0.08_{-0.12}^{+0.12})$.

Figure 4: Best-fit chemo-dynamical distributions for samples with different
vertical height cuts. We vary the vertical height constraints from $|z|>2.5$
to $4.0$ kpc, while the radial range is the same i.e. $r\in[7.5,8.5]$ kpc. The
columns from the left to the right are distributions of
$v_{r},v_{\theta},v_{\phi},\rm\rm{[Fe/H]}$, respectively. The rows from the
top to the bottom are for the disk, substructure and halo components,
respectively. The distributions are nearly invariant as $z_{\rm{cut}}$, except
for the $v_{\phi}$ distribution of the disk component, which shifts to a lower
velocity as the sample goes far away from the mid-plane.

###  3.2 Spatial Variation

We applied a vertical constraint of $|z|>2.5~{}$kpc to avoid the strong
dominance of the disk component in lower vertical region. In this Section, we
explore the spatial variation of our results by applying different radial and
vertical cuts to the samples to check if the results could be extrapolated to
the disk mid-plane.

####  3.2.1 Vertical Variation

To inspect the vertical variation, we select stars within $r\in[7.5,8.5]$ kpc,
and apply different vertical cuts of $|z|>z_{\rm{cut}}$ , varying
$z_{\rm{cut}}$ from 2.5 to 4 kpc. Fig. 4 displays the best-fit models for the
different vertical samples, with different rows representing the different
components and different columns for the chemo-dynamical distributions. Most
distributions of the three components do not vary much with the vertical cuts,
except for a decreasing azimuthal rotation with an increasing $z_{\rm{cut}}$,
which is expected. The radial velocity of the substructure shows slightly
stronger bimodal distribution as $z_{\rm{cut}}$ increases. Fig. 5 shows the
fractional contributions of the different components. The fraction of the disk
decreases with $z_{\rm{cut}}$ as expected, while the contributions of the
substructure and halo components both increase. Nevertheless, the relative
fraction of substructure component to the non-disk components are nearly
unchanged with $\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})\sim 0.85$.

Figure 5: Left: fractional contribution of different components as a function
of vertical height cut for samples within $r\in[7.5,8.5]~{}\rm{kpc}$ &
$|z|>z_{\rm cut}$. The green, blue and red dots are for the disk, substructure
and halo components, respectively. Right: fraction of the substructure
relative to the non-disk components. The error-bars represent the 2.1%, 50%
and 97.1% values of the MCMC results.

####  3.2.2 Radial Variation

Similarly, we explore the radial variation of the results by selecting
subsamples with the same vertical cut of $|z|>2.5$ kpc and different radial
cuts, varying from 7.5 to 10.0 kpc with a bin size of 0.5 kpc. The results are
shown in Figs. 6 and 7. As shown in Fig. 6, the chemo-dynamical distributions
of the disk and substructure are invariant for various radial ranges, while
the halo component shows some variation, especially in the metallicity. In the
range of $r\in[8.5,9.0]$ & $[9.5,10]$ kpc, the mean metallicity is much higher
with $\mu_{\rm{[Fe/H]}}\sim-1.5$ dex which is the upper limit of the prior and
also a larger dispersion than the other three radial bins. This metallicity
distribution is more similar to that of the substructure component, which
could be due to the strong degeneracy between the two components. In addition,
the halo fraction is always less than 10% in our samples, which means that its
chemo-dynamical distributions could be dominated and influenced by the other
two components. Nevertheless, as shown in the right panel of Fig. 7, the
relative fraction of substructure and halo still varies within 5% around
$\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})\sim 0.85$.

###  3.3 Parameter Degeneracy

There are 35 free parameters in total for the three components in the GMM. It
is important to inspect the degeneracy between the different parameters and
different components. The posterior distributions of these parameters except
for the fractions ($\mathcal{Q}_{h},\ \mathcal{Q}_{d},\ \mathcal{Q}_{s}$) are
illustrated in the Appendix, for each component separately. There is strong
degeneracy between the mean and dispersion of the metallicity for all three
components. More sepecifically, $\mu_{\rm{[Fe/H]}}$ and $\sigma_{\rm{[Fe/H]}}$
of the disk and substructure components show strong negative correlations,
while they show a positive correlation for the halo component. We also find
negative correlations between the mean and dispersion of the azimuthal
rotation of the disk component ($\mu_{\phi,d}$ vs. $\sigma_{\phi,d}$) and of
the radial motion of the substructure component ($\mu_{r,s}$ vs.
$\sigma_{r,s}$).

As the disk component is significantly different from the other two
components, the latter two could have a strong degeneracy, especially between
the radial motion and metallicity. Thus, we show the posterior distributions
of the parameters of $\mu_{r}$, $\sigma_{r}$, $\mu_{\rm{[Fe/H]}}$,
$\sigma_{\rm{[Fe/H]}}$ for the substructure and halo components, as well as
the fraction parameters in Fig. 8. Many correlations exist. (1) As expected
both the fractions of the disk ($\mathcal{Q}_{d}$) and the halo
($\mathcal{Q}_{h}$) components show negative correlations with the fraction of
the substructure ($\mathcal{Q}_{s}$) component, though these two show no
correlation. (2) The fraction of the halo component $\mathcal{Q}_{h}$ shows
strong correlations with the parameters of the metallicity of the substructure
and halo components. This is due mainly to that the halo component in our
sample contributes quite less with $\mathcal{Q}_{h}<0.1$. Its lower
metallicity leads it to emerge in the poor end of the metallicity distribution
in our sample. (3) The mean radial velocity of the substructure $\mu_{r,s}$
displays slight correlations with the fraction parameters $\mathcal{Q}_{d}$
and $\mathcal{Q}_{s}$. (4) The parameters of metallicity of the substructure
and halo components show strong correlations with each other. (5) The relative
fraction of the substructure in the non-disk components
$\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})$ is strongly effected by
the metallicity parameters. Though its value is not effected by the fraction
of disk component $\mathcal{Q}_{d}$, its determination is strongly influenced
by the poor end ([Fe/H] $<-1$) of the metallicity distribution of the sample.
It reveals the importance of the metallicity distribution of the sample for
separating these two components.

##  4 Angular Momentum Constraints

Figure 6: Same as Fig. 5 but for samples with the same vertical cut ($|z|>2.5$
kpc) and different radial ranges.

Closer to the mid-plane, the disk component is more dominant and the
identification of the substructure and halo components is more difficult. In
order to reduce the contribution and influence of the disk component, we apply
an additional angular momentum ($L_{z}$) constraint to our samples, and thus
attempt to utilize the GMM fitting to the region closer to the disk mid-plane.
As the substructure and halo components have slightly prograde rotation with
quite small mean rotation of $\mu_{\phi}\sim 20$ km s-1, we utilize the
angular momentum cut of $|L_{z}|<500$ kpc km s-1(Feuillet et al., 2020;
Belokurov et al., 2022), corresponding to $v_{\phi}\sim 62.5$ km s-1 at $R=8$
kpc.

Before moving closer to the mid-plane, we apply the angular momentum cut to
the sample within $r\in[7.5,8.5]~{}$kpc & $|z|>2.5~{}$kpc, and test the
performance of the GMM. The best-fit chemo-dynamical distributions of this
constrained sample are shown in Fig. 9. Fig. 10 shows the comparison of the
best-fit models with (dotted lines) and without (solid lines) the angular
momentum cut, with the 2.1%, 50% 97.9% percentiles of posterior distributions
of parameters in the GMM listed in Tables. 2-4 for the three components,
respectively.

Figure 7: Same as Fig. 5 but for samples with a same vertical cut ($|z|>2.5$
kpc) and different radial range.

The $v_{\phi}$ distribution of all three components have smaller
$\sigma_{\phi}$ as expected. The prograde rotation of the non-disk components
still remains, though the mean rotation $\mu_{\phi}$ decreases from $\sim 20$
km s-1 to $\sim 4$ km s-1 and $\sim 14$ km s-1 for the substructure and halo
components, respectively, as shown in Tables. 3 and 4. For all three
components, the medians of the best-fit distributions of $v_{r}$,
$v_{\theta}$, [Fe/H] are similar for the results with and without the angular
momentum constraint, as shown in Fig. 10 and Tables. 2-4. The angular momentum
constraint significantly decreases the contribution of the disk component and
the substructure becomes dominant, with component fractions of
$\mathcal{Q}_{d}=0.36^{+0.04}_{-0.03}$, $\mathcal{Q}_{s}=0.55^{+0.05}_{-0.05}$
and $\mathcal{Q}_{h}=0.09^{+0.04}_{-0.03}$. The resultant relative fraction of
the substructure and halo components is
$\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})=0.86^{+0.05}_{-0.06}$.
However, as the $v_{\phi}$ distributions of these two components have been
changed, the relative fractions of stars removed could be different. In order
to quantify and correct such a bias, we apply the angular momentum constraint
to the best-fit distribution obtained from MCMC, and thus deduce the relative
removed fractions in the substructure and halo components, respectively. After
correction, $\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})\sim 0.8$, which
is consistent with the results in Section 3.1.

Then we inspect the variation of the results with the vertical cut $z_{\rm
cut}$ for samples with the $L_{z}$ cut and additionally study two vertical
regions closer to the mid-plane. We select stars within $r\in[7.5,8.5]$ kpc &
$|L_{z}|<500$ kpc km s-1, and apply different vertical cuts of
$|z|>z_{\rm{cut}}$, varying $z_{\rm{cut}}$ from 0.5 to 4 kpc, including two
smaller $z_{\rm cut}$ of 0.5 and 1.5 kpc222With the angular momentum
constraint, the samples are almost the same for applying the vertical cuts
$z_{\rm{cut}}=0$ or $0.5$ kpc in this radial range. Therefore the case of
$z_{\rm{cut}}=0$ kpc is almost the same as $z_{\rm{cut}}=0.5$ kpc.. Fig. 11
shows the vertical variation of the best-fit chemo-dynamical distributions for
the $L_{z}$-cut samples similar to Fig. 4. The angular momentum constraint
cuts mainly the high rotation part of the disk component, which is more
significant at lower vertical height. The remaining stars in samples with
different $z_{\rm{cut}}$ show almost invariant chemo-dynamical distributions
for all three components, even for the two samples closer to the mid-plane. As
shown in Tabel. 4, for larger vertical cuts e.g. $z_{\rm{cut}}=3.5,4.0$ kpc,
the parameters of the halo component have larger uncertainties when the
$L_{z}$ constraint applied, which could be caused by the degeneracy between
the substructure and halo components. The radial velocity distribution of the
substructure component also shows slightly stronger bimodal distribution with
increasing $z_{\rm{cut}}$. As the velocity correlation $|\rho_{r\phi}|$ &
$|\rho_{\theta\phi}|\lesssim 0.1$, we can regard $v_{\phi}$ to be independent
on $v_{r}$ & $v_{\theta}$ for each component. Therefore, though the $v_{\phi}$
distribution has been significantly changed by the $L_{z}$ cut, the near
invariance of the results with $z_{\rm{cut}}$ means that we can safely
extrapolate our results to the Galactic plane.

Fig. 12 shows the comparison of fractions of the substructure component
relative to the non-disk components for results with (red) and without (blue)
the angular momentum constraint. The results from Necib et al. (2019a) are
also shown as the purple dots. This fraction is consistent for samples with
different vertical and radial cuts, and also for samples with or without the
angular momentum constraint. It demonstrates that for our sample
$\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})=0.8-0.85$, which is
slightly higher than that obtained in Necib et al. (2019a) of
$\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})=0.6-0.8$. Note that the
actual $\mathcal{Q}_{s}/(\mathcal{Q}_{s}+\mathcal{Q}_{h})$ ratio that Necib et
al. (2019a) adopted to deduce the DM velocity distribution is $76\%$, a value
computed using all stars with $d_{\odot}<4$ kpc & $|z|>2.5$ kpc.

##  5 Dark Matter Velocity Distribution and Detection Limits

Figure 8: Posterior distributions of possible degenerate model parameters from
the MCMC results of the sample in the region $r\in[7.5,8.5]$ kpc & $|z|>2.5$
kpc. The dashed lines of each 1D distribution represent the 2.1%,50%,97.9%
percentiles.

In this section we construct the local DM velocity distribution for the dual
DM model (i.e. the substructure and halo components), and estimate its
influence on the DM direct detection limits compared to the SHM. In Section
5.1 we present the conversion of the relative fraction of the substructure
component in the observed non-disk stellar components to the fraction of the
substructure component in the total DM, based on which we construct the local
DM velocity distribution. In Section 5.2 we show the impacts on the direct
detection limits where the spin-independent nuclear scattering is considered.

Figure 9: The best-fit chemo-dynamical distributions of the sample in the
region of $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc and $|L_{z}|<500$ kpc km s-1.

###  5.1 DM Velocity Distribution

In the SHM the DM halo is spherical and isotropic following the Maxwell-
Boltzmann distribution in the Galactocentric frame with the circular velocity
of the local standard of rest (LSR) $v_{0}=238~{}$km s-1 (Schönrich, 2012).
The velocity distribution is truncated at the escape velocity $v_{esc}=544$ km
s-1 at the Sun based on the work of Smith et al. (2009) to construct finite
extent halo which is also applied to the following dual halo model.

The dual halo consists of a DM substructure and a virialized SHM halo. In
order to construct the local DM velocity distribution from the stellar sample,
we follow the framework in Necib et al. (2019b) from which we can reproduce
their $c_{s}~{}\&~{}c_{h}$ in the following equation including the error
ranges. The velocity distribution of the dual halo model is:

$f_{\rm{dm}}(v)=\\!N\\!\left[\mathcal{Q}_{\ast,h}f_{h}(v)+\frac{c_{s}}{c_{h}}\mathcal{Q}_{\ast,s}f_{s}(v)\right]\
,$ (5)

where $N$ is the normalization factor,
$\mathcal{Q}_{\ast,j}=\frac{\mathcal{Q}_{j}}{\mathcal{Q}_{s}+\mathcal{Q}_{h}}\,(j=s,h)$,
and $f_{j}(v)$ is the best-fit velocity distribution given by the MCMC results
in Section 3.1, based on the correlation between the accreted stars and DM.
$\mathcal{Q}_{j}$ utilized in the following computation is also from Section.
3.1. $c_{j}=\frac{M_{\rm{peak}}}{M_{\ast,\rm{total}}}$, in which
$M_{\rm{peak}}$ and $M_{\ast,\rm{total}}$ are the peak halo mass and the total
stellar mass of the progenitor galaxy for each DM component, respectively.
However, $M_{\rm{peak}}$ and $M_{\ast,\rm{total}}$ of the progenitor galaxies
for the substructure and halo components are not given by observations
directly. Therefore we need to correlate $M_{\rm{peak}}$ and
$M_{\ast,\rm{total}}$ with observables like [Fe/H]. We utilize the SMHM
(stellar-mass-halo-mass) relation and the [Fe/H]-$M_{\ast,\rm{total}}$
(metallicity-stellar-mass) relation to derive the correlation between
$\left<\rm{[Fe/H]}\right>$ and $c_{j}$ as the following.

First, we generate $M_{\ast,\rm{total}}$ from $M_{\rm{peak}}$ by a Monte Carlo
procedure based on the SMHM relation. As illustrated in Behroozi et al.
(2013), this relation shows dependence on the redshift $z$, and different
trends for $M_{\rm peak}$ larger or smaller than a critical mass $M_{1}\sim
10^{11.5}M_{\odot}$. The general formalism is as below:

$\log(M_{\ast,\rm{total}})\\!=\\!\log_{10}(\epsilon
M_{1})\\!+\\!f\\!\\!\left[\log_{10}\\!\\!\left(\frac{M_{\rm
peak}}{M_{1}}\right)\\!\\!\right]\\!-\\!f(0)\ ,$ (6)

where

$f(x)=-\log_{10}(10^{-\alpha x}+1)+\delta\frac{(\log_{10}(1+\rm{exp}(\it
x)))^{\gamma}}{1+\rm{exp}(10^{-\it x})}\ .$ (7)

Figure 10: The comparison of best-fit chemo-dynamical distributions of the
sample in the region of $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc with (dotted
lines) and without (solid lines) the $|L_{z}|<500$ kpc km s-1 cut. Unlike
previous figures, the total area of each marginalized distribution (colored
curves) is normalized to the same constant in each panel.

For any given $M_{\rm peak}$, $M_{\ast,\rm{total}}$ is randomly selected from
the lognormal distribution with a median given by Eq. 6 and a scatter
$\sigma$. As illustrated in Garrison-Kimmel et al. (2016), for $M_{\rm
peak}>M_{1}$, free parameters in Eqs. 6 & 7, i.e.
$\epsilon,\alpha,\delta,\gamma$ and the scatter $\sigma$, are fixed to the
values given in Behroozi et al. (2013). While for $M_{\rm peak}\leq M_{1}$,
$\alpha\ \&\ \sigma$ are determined by the factor $\nu$ as:

$\alpha_{\nu}=0.25\nu^{2}-1.37\nu+1.69\ ,$ (8)
$\sigma_{\nu}=0.2+\nu\times(\log_{10}M_{\rm peak}-\log_{10}M_{1})\ .$ (9)

We set $\nu=-0.1$ in this benchmark model as in Necib et al. (2019b).

Secondly, we utilize the [Fe/H]-$M_{\ast,\rm{total}}$ relation for dwarf
galaxies in the MW investigated in Kirby et al. (2013) to generate
$\left<\rm{[Fe/H]}\right>$:

$\displaystyle\left<\rm{[Fe/H]}\right>$ $\displaystyle=(-1.69\pm 0.04)$ (10)
$\displaystyle+(0.30\pm
0.02)\log_{10}\frac{M_{\ast,\rm{total}}}{10^{6}~{}M_{\odot}}\ .$

$\left<\rm{[Fe/H]}\right>$ is randomly selected from the normal distribution
with a median given by Eq. 10 and a dispersion of 0.17.

Finally, we repeat the process above to compute the uncertainties of $c_{s}\
\&\ c_{h}$ by selecting points around $\left<\rm{[Fe/H]}\right>\sim-1.15$ and
$\left<\rm{[Fe/H]}\right>\sim-2.04$, representing the substructure and halo
components, respectively. This benchmark model gives
$\frac{c_{s}}{c_{h}}=0.05^{+0.10}_{-0.03}$ (16%, 50%, 84% percentiles). Thus
we can construct the Galactocentric velocity distribution of the local DM from
Eq. 5 as follows:

$f_{\rm{dm}}(v)=\mathcal{Q}_{\rm{dm,s}}f_{\rm{s}}(v)+\mathcal{Q}_{\rm{dm,h}}f_{\rm{h}}(v)\
,$ (11)

Figure 11: Same as Fig. 4 but for samples within $r\in[7.5,8.5]$ kpc and an
additional $|L_{z}|$ cut of $|L_{z}|<500$ kpc km s-1. The distributions are
roughly invariant even for the $v_{\phi}$ distribution of the disk component,
as the angular momentum constraints excludes most disk stars and reduces its
$v_{\phi}$ to a small value.

where $\mathcal{Q}_{\rm{dm,s}}$ and $\mathcal{Q}_{\rm{dm,h}}$ are the mass
fractions of DM in the substructure and isotropic halo components,
respectively. In Eq. 5 the normalization requires

$N(\mathcal{Q}_{\ast,h}+\frac{c_{s}}{c_{h}}\mathcal{Q}_{\ast,s})=1.$ (12)

From the above equation,
$\mathcal{Q}_{\rm{dm,s}}=N\times\frac{c_{s}}{c_{h}}\mathcal{Q}_{\ast,s}=25^{+24}_{-15}\%$
and $\mathcal{Q}_{\rm{dm,h}}=1-\mathcal{Q}_{\rm{dm,s}}$, where
$\mathcal{Q}_{\ast,s}=0.87$. The left panel of Fig. 13 shows the heliocentric
speed
($v_{\rm{helio}}=|\boldsymbol{v}_{\rm{Gal}}-\boldsymbol{v}_{\odot,\rm{Gal}}|$)
distribution of the substructure (blue dashed line) and halo (red dashed line)
weighted by their mass fractions, and their combination (black solid line).
For the comparison, the SHM is shown by the gray dashed line. Necib et al.
(2019a) shows that the heliocentric speed distribution of dual model
(consisting of the substructure and halo) shifts to the lower speed and has a
sharper peak compared with the SHM, shown as the dark-blue shaded region in
Fig. 13. Our result agrees with them but shifts to an even lower speed by
$\sim 7\%$. The uncertainty of the heliocentric speed distribution, shown as
the gray shaded region in Fig. 13, is contributed by two parts: one is the
uncertainty of the fractions of each component, the other is the uncertainty
in the velocity distribution especially for the halo component.

Figure 12: Left: the vertical variation of the fraction of the substructure
component relative to the non-disk components with $r\in[7.5,8.5]$ kpc. Right:
similar to the panel but for the radial variation.

###  5.2 Experimental Detection Limits

Considering the scattering between the DM particle and the target nuclear with
mass $m_{\chi}$ and $m_{T}$ respectively, the differential rate of scattering
events is given by the following equation:

$\frac{\rm{d}\it R}{\rm{d}\it
E_{R}}=\frac{\rho_{\chi}}{m_{\chi}m_{T}}\int_{v>v_{\rm{min}}}v\tilde{f}(\boldsymbol{v})\frac{\rm{d}\sigma}{\rm{d}\it
E_{R}}\rm{d}^{3}\it v\ ,$ (13)

where $E_{R}$ is the recoil energy deposited by the scattering, $\rho_{\chi}$
is the local DM density set as $0.3~{}\rm Gev/cm^{3}$ (Amsler et al., 2008)
for both the SHM and the dual model, and $\tilde{f}(\boldsymbol{v})$ is the
heliocentric velocity distribution333Critically, it should be the laboratory
velocity distribution considering cyclically velocity of the earth with
respect to the Sun. To simplify the computation, the time average of
laboratory velocity over one period is utilized here, namely the heliocentric
velocity. of the local DM. $v_{\rm{min}}$ is the minimum velocity for a DM
particle depositing the recoil energy $E_{R}$ defined as:

$v_{\rm{min}}=\sqrt{\frac{m_{T}E_{R}}{2\mu^{2}}}$ (14)

where $\mu$ is the reduced mass of the DM particle and target nuclear. The
differential scattering cross section $\frac{\rm{d}\sigma}{\rm{d}\it E_{R}}$
is defined as:

$\frac{\rm{d}\it\sigma}{\rm{d}\it
E_{R}}=\frac{m_{T}A^{2}\sigma_{\rm{SI}}}{2\mu^{2}v^{2}}F^{2}(E_{R})\ ,$ (15)

where $A$ is the atomic number of the target nuclear, $\sigma_{\rm{}_{SI}}$ is
the spin-independent scattering cross section and $F(E_{R})$ is the nuclear
form factor set as the Helm form (Helm, 1956). Then Eq. 13 can be written as:

$\frac{\rm{d}\it R}{\rm{d}\it
E_{R}}=\frac{\rho_{\chi}A^{2}\sigma_{\rm{SI}}}{2m_{\chi}\mu^{2}}F^{2}(E_{R})g(v_{\rm
min})\ ,$ (16)

where $g(v_{\rm min})$ is defined as:

$g(v_{\rm{min}})=\int_{v>v_{\rm{min}}}\frac{\tilde{f}(\boldsymbol{v})}{v}\rm{d}^{3}\it
v\ .$ (17)

Finally the total number of events is given by:

$N_{e}=MT\int_{0}^{\infty}\phi(E_{R})\frac{\rm{d}\it R}{\rm{d}\it
E_{R}}\rm{d}\it E_{R}\ ,$ (18)

where $M$ is the total mass of the target, $T$ is the exposure time and
$\phi(E_{R})$ is the detector acceptance, i.e. the efficiency for different
recoil energy, which differs for different experiments.

For a specific experiment, the number of events is the function of $m_{\chi}$,
$\sigma_{\rm SI}$ and $\tilde{f}(\boldsymbol{v})$. As an example, we utilize
corresponding parameters of the Xenon direct detection experiment PandaX-4T
(Meng et al., 2021) with an exposure of 0.63 tonne $\times$ year. Applying the
SHM, we convert the published exclusion curve of PandaX-4T into $N_{e}$ vs
$m_{\chi}$ with the efficiency curve given in Meng et al. (2021). Then
applying the DM velocity distribution deduced in Section. 5.1, we re-compute
the cross section limit resulting in the same $N_{e}$ vs $m_{\chi}$. All
points have $m_{\chi}\gtrsim 5$ Gev due to the low signal efficiency for low
mass DM particles.

The right panel of Fig. 13 shows the 90% C.L. direct detection limits of the
DM substructure (blue dashed), the DM halo (red dashed) and their total
contribution (black solid). When computing the $g(v_{\rm{min}})$ we use the
fraction weighted velocity distribution of the substructure and halo
components. Compared with the SHM (gray dashed), the dual model results in a
weaker limit for lower mass DM due to the softer DM kinetic energy
distribution, and a slightly stronger limit for DM mass above 150 GeV for
which the signal acceptance increases for a given set of selection cuts.

##  6 Discussions

Figure 13: Left: the heliocentric speed distributions of the best-fit
substructure (blue dashed line) and halo (red dashed line) components,
weighted by their relative fractions in the Solar neighbourhood. The
combination of the two is shown as the black solid line, while the gray shaded
region displays the uncertainty. The gray dashed line represents the SHM. The
dark-blue shaded region shows the result from Necib et al. (2019a). Right: 90%
C.L. upper DM direct detection limits on the cross section ($\sigma_{\rm SI}$)
vs. DM mass ($m_{\chi}$) plane, calculated based on the DM direct detection
experiment PandaX-4T (Meng et al., 2021). The representation of each line is
the same as the left panel.

###  6.1 Methods of Decreasing the Contribution of Disk Stars

In this work our main interest is the non-disk stars, i.e. the substructure
and main halo components. The disk-dominance could lead to uncertainty for
other two components. There are many methods to remove the disk stars, such as
constraints on the spatial volume, the chemical abundance, the kinematic
information etc. Here we try to decrease the contribution of disk stars
through different methods and assess the influence.

Wu et al. (2022) applied an effective selection criteria ($|z|>5$ kpc &
${\rm{[Fe/H]}}<-0.5$) to remove possible disk stars. As our interested region
is the Solar neighbourhood with lower vertical height, we attempt to cut the
stellar sample with a more strict criteria of [Fe/H] $\lesssim-1$ to exclude
most of possible disk stars. We then fit the sample using a two-component
model with only the substructure and halo components. However, as the mean
[Fe/H] of the substructure component is $\sim-1.2$ dex, this metallicity cut
would also exclude nearly half of the substructure stars. The GMM method
failed to fit the stellar sample, as the peak metallicity of the substructure
component would prefer a value beyond the upper metallicity boundary of $-1$
dex.

In Section 4 we show one of the possible kinematic constraints to decrease the
contribution of disk, which supports our results in Section 3 and also
extrapolates our results to the Galactic plane. Another commonly applied
kinematic cut to select the GES stars is the Toomre diagram as shown in Fig.
14. After excluding disk stars from the Toomre diagram, the chemo-dynamical
distribution of our sample shows a non-negligible contribution of metal-rich
stars with peak [Fe/H] $\sim-0.6$ dex. Besides, $v_{r}\ \&\ v_{\phi}$
distributions of the remained sample show prominent asymmetry as the cut is
asymmetric in velocity space, which excludes more stars with prograde
rotation. Thus, the GMM again failed to fit the sample after the Toomre
diagram cut due to such an asymmetry. Though we can separate the bimodal
Gaussian distribution into two Gaussian distributions with different
fractions, the fitting only gives consistent results with Section 3 in the
means of velocities. It may bias the resultant velocity distributions of the
non-disk components.

The asymmetry of $v_{r}\ \&\ v_{\phi}$ distribution indicates that the meal-
rich peak might correspond to the hot thick disk and the Nyx stream (Necib et
al., 2020). Furthermore, Donlon & Newberg (2023) combined this kinematic cut
with an additional chemical abundance cut (Das et al., 2020; Buder et al.,
2021; Belokurov & Kravtsov, 2022) and obtained a local stellar halo sample.
Their selection requires more chemical abundance information such as
$[\rm{Al/Fe}],\ [\rm{Mg/Mn}],\ [\rm{C/Fe}]$ besides [Fe/H]. In the future,
based on more chemical abundance information from LAMOST, we might obtain more
detailed structures of the local stellar halo.

Figure 14: Toomre diagram for the sample in the region of $r\in[7.5,8.5]$ kpc
$\&$ $|z|>2.5$ kpc. The black circle represents
$|\boldsymbol{v}-\boldsymbol{v}_{lsr}|=210$ km s-1 where
$\boldsymbol{v}_{lsr}=(0,238,0)$ km s-1, which presumably separates halo stars
from disk stars.

###  6.2 Degeneracy Between Substructure and Halo

In Section 3.3 we investigate the potential degeneracy between the
substructure and halo components especially for the radial velocity and
metallicity distribution. As the halo component contributes a small fraction
to the total sample, we also tried a two-component model including only the
disk and substructure components to fit the sample, and compare with the
results from the three-component model.

Besides in the regions of $r\in[8.5,9.0]$ kpc and $r\in[9.5,10.0]$ kpc when
$|z|>2.5$ kpc, the halo component has a mean metallicity
$\mu_{\rm{[Fe/H]}}\sim-1.5$ dex. The different $\mu_{\rm{[Fe/H]}}$ from that
in other radial bins could be caused by the strong degeneracy with the
substructure component. Therefore we firstly apply the two-component model in
these two regions and compare with the tree-component model.

The Bayesian information criterion (BIC, Schwarz, 1978) is utilized to
quantify which model performs better. The BIC is defined as:

$\rm{BIC}=\it k\ln N-\rm 2\times\ln\mathcal{L}_{MAP}\ ,$ (19)

where $k$ is the number of free parameters in the model, $N$ is the size of
the sample and $\mathcal{L}_{\rm MAP}$ is the maximum of the likelihood
function deciding the best-fit parameters. Fitting to the same sample, these
two models share the same $N$ while the two-component model contains 23 free
parameters in total.

We calculate the BIC for each model and find
$|\rm{BIC_{2com}}-\rm{BIC_{3com}}|\lesssim 10$, which means both models fit
the sample statistically similarly. Nevertheless, for other radial bins,
$\rm{BIC_{2com}}-\rm{BIC_{3com}}\gtrsim 100$, which means that the three-
component model performs better and the contribution from the halo component
is not trivial.

One of possible contributions to the degeneracy is from other substructures
with similar chemical abundance and kinematic properties as the substructure
and halo components. As mentioned in Wu et al. (2022) the most significant
contaminant is the Sgr stream with $\mu_{\rm{[Fe/H]}}\sim-1.3$ dex (Yang et
al., 2019), which has similar metallicity as the substructure component.
Nevertheless, as we focus on the Solar neighbourhood with $r<10$ kpc, the
contribution of the Sgr stream is believed to be negligible (Yang et al.,
2019). Their results with the Sgr stream removed showed that the radial motion
of the GES component, i.e. the substructure component in our work, is affected
by other substructures. With most substructures removed, the radial velocity
distribution of the GES component shows stronger bimodal distribution. Note
that, the hot thick disk and the Nyx stream discussed in Necib et al. (2020)
might contribute to the disk and the substructure components, which need more
careful studies in the future.

###  6.3 Limitation of Tracing DM Velocity Distribution with Stars

The correlation in velocity distribution between the accreted stars and DM is
the basic assumption in this work, which is still being debated. Such a
correlation is found for those accreted from luminous satellite galaxies.
However, Wang et al. (2011) and Necib et al. (2019a) showed that a non-
negligible fraction of DM of Milky Way-like galaxies could be accreted from
smooth accretion of dark satellites, which cannot be traced by accreted stars.
This limitation is also mentioned in Bozorgnia et al. (2019) to explain the
weak correlation in the velocity distributions of metal-poor stars and DM
particles in Auriga simulations. We tried to investigate such correlation in
TNG-50 simulations (Pillepich et al., 2019; Nelson et al., 2019) and found a
rough correlation between stars and DM accreted from the same satellite.
However, the resolution of TNG-50 simulations is still not enough for most
Milky Way-like galaxies to check the velocity distribution for accreted
particles in the Solar neighbourhood. Higher resolution for Milky Way-like
galaxies is needed to further investigate the correlation between the stars
and DM. This correlation may also depend on the merger orbit of the satellite
galaxy.

The benchmark model of converting the relative contribution of the
substructure in the stellar sample to that in DM is very sensitive to the
metallicity of each stellar component. Different spectroscopic surveys have
different selection effects thus different metallicity distribution of the
stellar sample. It may contribute to the difference between our results and
Necib et al. (2019a). Besides, this conversion also depends on the SMHM
relation of the satellites which varies for different simulations. All these
could result in additional uncertainties in the relative contribution from the
DM substructure.

Besides the contribution of the radial merger event of GES happened around
8$-$10 Gyr ago (Sahlholdt et al., 2019; Bignone et al., 2019; Bonaca et al.,
2020), other recent merger events may also have non-negligible influence on
the local DM density and velocity distribution. O’Hare et al. (2020) computed
the impact of the S1 stream on WIMP and axion detectors, which passes through
the Solar neighbourhood on a low inclination and counter-rotating orbit, with
a progenitor mass comparable to the present-day Fornax dwarf spheroidal.
Donlon et al. (2022) identified three radial merger events, i.e. Virgo radial
merger, Nereus and Cronus, using local halo dwarf stars. They thought stars
comprising the GES velocity structure are a combination of those components.
Smith-Orlik et al. (2023) studied the effect of the Large Magellanic Cloud
(LMC) on the DM distribution in the Solar neighbourhood utilizing simulations
of Milky Way analogues. They found that DM particles in the Solar
neighbourhood originating from the LMC analogue dominate the high velocity
tail of the local DM velocity distribution, and the native DM particles of the
MW in the Solar vicinity are boosted to higher speeds as a result of a
response to the LMC’s motion.

##  7 Conclusion

In this work, we attempt to study the local DM velocity distribution based on
the correlations in velocity distribution between stars and DM belonging to
different components. We apply the GMM to separate the disk, the halo
substructure and the main halo components in the chemo-dynamical space
utilizing K giants of LAMOST DR8 cross-matched with Gaia DR3. The substructure
component is radially anisotropic possibly related to the GES debris flow,
while the halo component is referred to as the isotropic halo accreted
earliest and thus virialized. Our main results are derived from the subsample
within $r\in[7.5,10]$ kpc & $|z|>2.5$ kpc. We also explore the spatial
variation of the best-fit distributions. We try to apply an additional angular
momentum constraint to reduce the dominance of the disk component. The main
results are consistent and are invariant with the vertical cut, which verifies
the robustness of extrapolating our main results to the Solar vicinity. Based
on the stellar chemo-dynamical distributions, we compute the velocity
distribution of the dual DM, i.e. the substructure and halo components,
applying the [Fe/H]-$M_{\ast,\rm{total}}$ relation and the SMHM relation. This
modified DM velocity distribution is different from the SHM, resulting in
different exclusion curves ($\sigma_{\rm SI}$ vs. $m_{\chi}$) for DM direct
detection experiments.

Our main results are as follows:

  * •

The GMM performs well in separating stars belonging to different components in
the chemo-dynamical space. For our K giant stars within $r\in[7.5,8.5]$ kpc
$\&$ $|z|>2.5$ kpc, the substructure component, i.e. the GES debris flow,
contributes about 87% of the non-disk stars
($\mathcal{Q}_{s}/(\mathcal{Q}_{h}+\mathcal{Q}_{s})=0.87^{+0.03}_{-0.04}$).
This dominant stellar halo component is highly radially anisotropic with an
anisotropy parameter
$\beta=1-\frac{\sigma_{\theta}^{2}+\sigma_{\phi}^{2}}{2\sigma_{r}^{2}}=0.71$,
which is smaller than some previous works (e.g. Belokurov et al., 2018; Myeong
et al., 2018; Bird et al., 2019; Wu et al., 2022). The halo component referred
to as the virialized halo is more isotropic with $\beta=0.23$. Both components
show a slightly prograde rotation with $\mu_{\phi}\sim 20$ km s-1. The median
metallicities for the substructure and halo components in our sample are
$-1.15$ and $-2.04$, respectively. The substructure component is slightly more
metal-rich and the halo component is more metal-poor than the corresponding
compnents in Necib et al. (2019a) and Wu et al. (2022).

  * •

The best-fit chemo-dynamical distributions of the non-disk components are
almost invariant for subsamples with different vertical cuts or different
radial ranges in $r\in[7.5,10]$ kpc and $|z|>2.5$ kpc. Furthermore we
extrapolate our results to the Galactic plane with an additional angular
momentum ($L_{z}$) constraint to reduce the dominance of the disk component.
The main results with the angular momentum constraint, e.g.
$\mathcal{Q}_{s}/(\mathcal{Q}_{h}+\mathcal{Q}_{s})$, are consistent with those
from samples without the $L_{z}$ cut. These tests demonstrates the robustness
of extrapolating our main results to the Solar vicinity.

  * •

Some parameters in the GMM are degenerate, especially between the parameters
of metallicity ($\mu_{\rm{[Fe/H]}}$, $\sigma_{\rm{[Fe/H]}}$) and the fractions
($\mathcal{Q}_{s}$, $\mathcal{Q}_{h}$) of the substructure and halo
components. These correlations imply the degeneracy between the substructure
and halo components, which is partly due to the relatively small contribution
of the halo component. Nevertheless, the BIC test indicates that in most
regions (except r$\in[8.5,9.0]$ kpc and $[9.5,10]$ kpc), the three-component
model is preferred over a two-component model, which shows that the halo
component is non-negligible.

  * •

Utilizing the [Fe/H]-$M_{\ast,\rm{total}}$ relation and the SMHM relation, we
convert the relative contribution of the substructure component to the total
halo stars ($\mathcal{Q}_{s}/(\mathcal{Q}_{h}+\mathcal{Q}_{s})\sim 85\%$) to
the relative contribution in DM of $25_{-15}^{+24}\%$ in the Solar
neighbourhood. As a comparison, Necib et al. (2019a) indicated the
substructure contributes 60%-80% of the local total halo stars and
$42^{+26}_{-22}\%$ of DM.

  * •

Considering the contributions from both the substructure and halo DM, the
local DM heliocentric velocity distribution shifts to the lower speed and has
a sharper peak compared to the SHM. For the Xenon direct detection experiment
PandaX-4T, the dual DM model indicates a significant larger cross section
$\sigma_{\rm{SI}}$ for low mass DM of $m_{\chi}\lesssim 150$ Gev.

## Acknowledgement

The research presented here is partially supported by the National Key R&D
Program of China under grant No. 2018YFA0404501; by the National Natural
Science Foundation of China under grant Nos. 11773052, 11761131016, 11333003,
12025302, 12103031; and by the “111” Project of the Ministry of Education
under grant No. B20019. This project was developed in part at the LAMOST-Gaia
Sprint 2022, supported by the National Natural Science Foundation of China
(NSFC) under grants 11873034 and U2031202. R.G. is supported by Initiative
Postdocs Supporting Program (No. BX2021183), funded by China Postdoctoral
Science Foundation. J.S. acknowledges support from a Newton Advanced
Fellowship awarded by the Royal Society and the Newton Fund. J.L. is supported
by NSFC (No. 1209060), the Hongwen Foundation in Hong Kong, and Tencent and
New Cornerstone Science Foundation in China. X.-X.X. and L.Z. acknowledge the
support from National Key Research and Development Program of China No.
2019YFA0405504, National Natural Science Foundation of China (NSFC) under
grants No. 11988101, CAS Project for Young Scientists in Basic Research grant
No. YSBR-062 and YSBR-092, and China Manned Space Project with No. CMS-
CSST-2021-B03.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectro-scopic
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 of China, 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.

This work made use of the Gravity Supercomputer at the Department of
Astronomy, Shanghai Jiao Tong University.

We provide more details of the GMM results in this appendix. In Tables. 2-4 we
provide posterior distribution of parameters of the GMM for the disk,
substructure and halo component, respectively. The conner plots for the sample
in the region $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc are displayed in Figs.
15-17.

Table 2: Parameters for the GMM of the disk component with and without $L_{z}$ cut in the region of $r\in[7.5,8.5]$ kpc & $|z|>z_{cut}$. $z_{cut}$ | $\mu_{r}$ | $\mu_{\theta}$ | $\mu_{\phi}$ | $\sigma_{r}$ | $\sigma_{\theta}$ | $\sigma_{\phi}$ | $\rho_{r\theta}$ | $\rho_{r\phi}$ | $\rho_{\theta\phi}$ | $\mu_{\rm{[Fe/H]}}$ | $\sigma_{\rm{[Fe/H]}}$
---|---|---|---|---|---|---|---|---|---|---|---
kpc | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 |  |  |  | dex | dex
With Angular Momentum Constraints
$0.5$ | $-2.42_{-5.52}^{+5.58}$ | $-2.58_{-3.57}^{+3.48}$ | $24.79_{-2.49}^{+2.32}$ | $94.72_{-4.72}^{+4.85}$ | $60.74_{-2.53}^{+2.62}$ | $34.11_{-1.87}^{+1.77}$ | $0.09_{-0.06}^{+0.06}$ | $-0.06_{-0.06}^{+0.06}$ | $-0.01_{-0.06}^{+0.06}$ | $-0.59_{-0.01}^{+0.01}$ | $0.16_{-0.01}^{+0.01}$
$1.5$ | $-1.38_{-6.76}^{+6.85}$ | $-3.33_{-4.04}^{+3.88}$ | $23.23_{-2.78}^{+3.03}$ | $95.13_{-5.75}^{+5.39}$ | $60.03_{-2.93}^{+2.92}$ | $35.72_{-2.15}^{+1.99}$ | $0.09_{-0.08}^{+0.07}$ | $0.04_{-0.07}^{+0.07}$ | $0.01_{-0.08}^{+0.07}$ | $-0.59_{-0.01}^{+0.01}$ | $0.15_{-0.01}^{+0.01}$
$2.5$ | $-4.97_{-7.75}^{+7.80}$ | $-3.41_{-4.33}^{+4.20}$ | $21.58_{-3.09}^{+3.16}$ | $97.84_{-6.46}^{+6.60}$ | $56.58_{-3.17}^{+3.31}$ | $38.09_{-2.35}^{+2.33}$ | $0.18_{-0.08}^{+0.08}$ | $0.05_{-0.08}^{+0.08}$ | $-0.01_{-0.08}^{+0.08}$ | $-0.60_{-0.02}^{+0.02}$ | $0.14_{-0.02}^{+0.02}$
$3.0$ | $-4.51_{-8.64}^{+8.42}$ | $-4.21_{-4.85}^{+4.89}$ | $20.82_{-3.50}^{+3.68}$ | $95.67_{-7.29}^{+7.57}$ | $55.17_{-3.61}^{+3.87}$ | $39.33_{-2.60}^{+2.64}$ | $0.19_{-0.09}^{+0.09}$ | $0.03_{-0.09}^{+0.09}$ | $-0.00_{-0.10}^{+0.09}$ | $-0.61_{-0.02}^{+0.02}$ | $0.14_{-0.02}^{+0.02}$
$3.5$ | $-5.62_{-9.79}^{+9.74}$ | $-4.74_{-5.24}^{+5.75}$ | $18.35_{-4.04}^{+4.42}$ | $95.76_{-7.97}^{+8.52}$ | $53.80_{-4.15}^{+4.41}$ | $41.68_{-3.08}^{+3.17}$ | $0.17_{-0.11}^{+0.11}$ | $0.01_{-0.11}^{+0.10}$ | $0.01_{-0.10}^{+0.11}$ | $-0.62_{-0.02}^{+0.02}$ | $0.14_{-0.02}^{+0.02}$
$4.0$ | $-1.29_{-11.86}^{+11.24}$ | $-6.58_{-6.40}^{+6.26}$ | $16.58_{-4.80}^{+5.22}$ | $96.07_{-9.38}^{+9.32}$ | $51.68_{-4.31}^{+4.82}$ | $44.03_{-3.64}^{+3.83}$ | $0.19_{-0.14}^{+0.13}$ | $0.05_{-0.14}^{+0.12}$ | $0.00_{-0.12}^{+0.11}$ | $-0.62_{-0.03}^{+0.03}$ | $0.15_{-0.03}^{+0.03}$
Without Angular Momentum Constraints
$2.5$ | $6.02_{-2.90}^{+2.78}$ | $1.37_{-1.73}^{+1.69}$ | $134.73_{-2.87}^{+2.74}$ | $81.24_{-2.33}^{+2.39}$ | $48.79_{-1.30}^{+1.34}$ | $65.89_{-2.17}^{+2.12}$ | $0.16_{-0.04}^{+0.04}$ | $-0.09_{-0.04}^{+0.04}$ | $0.02_{-0.04}^{+0.04}$ | $-0.55_{-0.01}^{+0.01}$ | $0.17_{-0.01}^{+0.01}$
$3.0$ | $4.10_{-4.03}^{+3.92}$ | $0.79_{-2.33}^{+2.29}$ | $114.24_{-3.61}^{+3.71}$ | $85.80_{-3.16}^{+3.20}$ | $50.45_{-1.82}^{+1.83}$ | $71.63_{-2.65}^{+2.74}$ | $0.21_{-0.05}^{+0.05}$ | $-0.04_{-0.05}^{+0.05}$ | $0.05_{-0.05}^{+0.05}$ | $-0.57_{-0.01}^{+0.01}$ | $0.16_{-0.01}^{+0.01}$
$3.5$ | $1.58_{-5.46}^{+5.70}$ | $-1.17_{-3.04}^{+3.21}$ | $95.67_{-5.20}^{+5.08}$ | $90.80_{-4.42}^{+4.15}$ | $52.65_{-2.48}^{+2.68}$ | $77.05_{-3.74}^{+3.73}$ | $0.27_{-0.06}^{+0.06}$ | $-0.05_{-0.07}^{+0.07}$ | $0.04_{-0.07}^{+0.06}$ | $-0.59_{-0.02}^{+0.01}$ | $0.16_{-0.01}^{+0.01}$
$4.0$ | $2.67_{-7.37}^{+7.61}$ | $-3.38_{-4.36}^{+4.25}$ | $74.88_{-6.60}^{+6.55}$ | $93.99_{-6.05}^{+6.30}$ | $54.09_{-3.42}^{+3.53}$ | $78.41_{-5.05}^{+5.16}$ | $0.28_{-0.08}^{+0.08}$ | $-0.08_{-0.08}^{+0.09}$ | $0.03_{-0.08}^{+0.08}$ | $-0.63_{-0.03}^{+0.03}$ | $0.18_{-0.03}^{+0.03}$
Table 3: Same as Table. 2 but for the substructure component. $z_{cut}$ | $\mu_{r}$ | $\mu_{\theta}$ | $\mu_{\phi}$ | $\sigma_{r}$ | $\sigma_{\theta}$ | $\sigma_{\phi}$ | $\rho_{r\theta}$ | $\rho_{r\phi}$ | $\rho_{\theta\phi}$ | $\mu_{\rm{[Fe/H]}}$ | $\sigma_{\rm{[Fe/H]}}$
---|---|---|---|---|---|---|---|---|---|---|---
kpc | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 |  |  |  | dex | dex
With Angular Momentum Constraints
$0.5$ | $132.92_{-7.85}^{+7.78}$ | $-5.34_{-3.92}^{+3.95}$ | $3.93_{-2.23}^{+2.28}$ | $110.25_{-5.10}^{+5.70}$ | $78.50_{-2.93}^{+3.01}$ | $42.41_{-1.49}^{+1.48}$ | $0.10_{-0.06}^{+0.07}$ | $0.05_{-0.06}^{+0.06}$ | $-0.02_{-0.05}^{+0.05}$ | $-1.22_{-0.04}^{+0.04}$ | $0.38_{-0.03}^{+0.03}$
$1.5$ | $137.26_{-7.91}^{+8.21}$ | $-4.59_{-4.26}^{+4.16}$ | $3.89_{-2.34}^{+2.26}$ | $109.69_{-5.20}^{+5.98}$ | $76.73_{-3.22}^{+3.07}$ | $42.96_{-1.56}^{+1.63}$ | $0.13_{-0.07}^{+0.07}$ | $-0.04_{-0.07}^{+0.06}$ | $0.01_{-0.06}^{+0.06}$ | $-1.23_{-0.04}^{+0.04}$ | $0.38_{-0.04}^{+0.04}$
$2.5$ | $141.71_{-8.44}^{+8.54}$ | $-6.08_{-4.64}^{+4.41}$ | $3.62_{-2.76}^{+2.79}$ | $107.35_{-5.54}^{+6.38}$ | $76.50_{-3.54}^{+3.78}$ | $44.59_{-1.87}^{+1.98}$ | $0.18_{-0.07}^{+0.07}$ | $-0.05_{-0.08}^{+0.08}$ | $0.03_{-0.06}^{+0.06}$ | $-1.25_{-0.04}^{+0.04}$ | $0.37_{-0.04}^{+0.04}$
$3.0$ | $144.97_{-9.20}^{+9.13}$ | $-8.35_{-5.14}^{+4.97}$ | $3.64_{-3.11}^{+2.94}$ | $105.25_{-5.59}^{+6.13}$ | $77.34_{-3.91}^{+4.21}$ | $45.75_{-2.21}^{+2.26}$ | $0.20_{-0.08}^{+0.08}$ | $-0.06_{-0.09}^{+0.09}$ | $0.03_{-0.07}^{+0.07}$ | $-1.28_{-0.06}^{+0.05}$ | $0.39_{-0.05}^{+0.05}$
$3.5$ | $145.87_{-10.52}^{+10.52}$ | $-9.78_{-5.61}^{+5.73}$ | $3.13_{-3.56}^{+3.40}$ | $105.07_{-6.10}^{+7.22}$ | $78.23_{-5.43}^{+5.59}$ | $46.93_{-2.57}^{+2.62}$ | $0.20_{-0.08}^{+0.08}$ | $-0.01_{-0.09}^{+0.09}$ | $0.02_{-0.08}^{+0.08}$ | $-1.32_{-0.07}^{+0.07}$ | $0.42_{-0.06}^{+0.04}$
$4.0$ | $148.74_{-11.53}^{+14.90}$ | $-10.77_{-6.90}^{+6.89}$ | $4.02_{-4.52}^{+3.83}$ | $106.73_{-7.79}^{+9.01}$ | $79.52_{-7.01}^{+5.82}$ | $48.81_{-3.05}^{+2.74}$ | $0.21_{-0.10}^{+0.09}$ | $0.00_{-0.10}^{+0.10}$ | $0.03_{-0.09}^{+0.09}$ | $-1.33_{-0.06}^{+0.07}$ | $0.43_{-0.06}^{+0.04}$
Without Angular Momentum Constraints
$2.5$ | $120.09_{-6.61}^{+6.44}$ | $-6.88_{-3.34}^{+3.35}$ | $22.59_{-4.61}^{+4.74}$ | $108.52_{-4.53}^{+5.20}$ | $79.88_{-2.67}^{+2.64}$ | $93.43_{-3.20}^{+3.20}$ | $0.15_{-0.05}^{+0.05}$ | $-0.08_{-0.06}^{+0.06}$ | $0.06_{-0.04}^{+0.04}$ | $-1.15_{-0.04}^{+0.04}$ | $0.42_{-0.02}^{+0.02}$
$3.0$ | $131.37_{-7.21}^{+6.73}$ | $-9.23_{-4.15}^{+3.84}$ | $15.73_{-5.31}^{+5.32}$ | $106.98_{-4.78}^{+5.31}$ | $81.67_{-3.05}^{+3.07}$ | $93.77_{-3.56}^{+3.56}$ | $0.18_{-0.06}^{+0.06}$ | $-0.09_{-0.06}^{+0.07}$ | $0.08_{-0.05}^{+0.05}$ | $-1.21_{-0.04}^{+0.04}$ | $0.41_{-0.03}^{+0.02}$
$3.5$ | $134.31_{-7.87}^{+7.57}$ | $-10.26_{-4.67}^{+4.93}$ | $13.73_{-5.80}^{+5.52}$ | $107.36_{-5.27}^{+6.02}$ | $83.01_{-3.77}^{+3.89}$ | $94.03_{-4.27}^{+3.95}$ | $0.17_{-0.07}^{+0.07}$ | $-0.10_{-0.07}^{+0.08}$ | $0.08_{-0.05}^{+0.06}$ | $-1.24_{-0.05}^{+0.04}$ | $0.40_{-0.03}^{+0.03}$
$4.0$ | $144.29_{-9.84}^{+10.57}$ | $-12.93_{-5.64}^{+5.53}$ | $9.98_{-7.28}^{+7.19}$ | $107.00_{-6.14}^{+7.01}$ | $86.96_{-4.28}^{+4.58}$ | $94.47_{-4.82}^{+4.89}$ | $0.18_{-0.08}^{+0.07}$ | $-0.11_{-0.09}^{+0.08}$ | $0.10_{-0.07}^{+0.07}$ | $-1.28_{-0.05}^{+0.05}$ | $0.39_{-0.04}^{+0.04}$
Table 4: Same as Table. 2 but for the halo component. $z_{cut}$ | $\mu_{r}$ | $\mu_{\theta}$ | $\mu_{\phi}$ | $\sigma_{r}$ | $\sigma_{\theta}$ | $\sigma_{\phi}$ | $\rho_{r\theta}$ | $\rho_{r\phi}$ | $\rho_{\theta\phi}$ | $\mu_{\rm{[Fe/H]}}$ | $\sigma_{\rm{[Fe/H]}}$
---|---|---|---|---|---|---|---|---|---|---|---
kpc | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 | km s-1 |  |  |  | dex | dex
With Angular Momentum Constraints
$0.5$ | $-10.30_{-18.62}^{+19.01}$ | $-0.34_{-15.32}^{+15.60}$ | $12.94_{-7.50}^{+7.15}$ | $119.41_{-15.62}^{+15.22}$ | $117.03_{-11.94}^{+15.41}$ | $51.20_{-5.01}^{+6.43}$ | $-0.04_{-0.14}^{+0.14}$ | $-0.10_{-0.16}^{+0.15}$ | $0.08_{-0.14}^{+0.13}$ | $-2.03_{-0.06}^{+0.06}$ | $0.21_{-0.04}^{+0.04}$
$1.5$ | $-11.23_{-21.51}^{+24.27}$ | $0.36_{-17.95}^{+17.47}$ | $14.23_{-8.38}^{+10.95}$ | $119.92_{-20.57}^{+18.02}$ | $115.16_{-11.84}^{+17.45}$ | $53.88_{-6.08}^{+8.08}$ | $-0.06_{-0.17}^{+0.15}$ | $0.10_{-0.18}^{+0.19}$ | $-0.04_{-0.15}^{+0.15}$ | $-2.03_{-0.06}^{+0.08}$ | $0.21_{-0.04}^{+0.06}$
$2.5$ | $-10.01_{-22.68}^{+21.43}$ | $7.05_{-17.90}^{+18.42}$ | $14.43_{-9.00}^{+10.15}$ | $123.97_{-17.48}^{+17.27}$ | $115.72_{-13.28}^{+18.78}$ | $55.50_{-6.12}^{+8.21}$ | $-0.05_{-0.16}^{+0.16}$ | $0.11_{-0.17}^{+0.17}$ | $-0.07_{-0.16}^{+0.16}$ | $-2.03_{-0.07}^{+0.08}$ | $0.22_{-0.05}^{+0.05}$
$3.0$ | $-11.83_{-32.01}^{+30.17}$ | $6.98_{-24.40}^{+25.84}$ | $16.08_{-11.39}^{+16.55}$ | $125.16_{-22.83}^{+21.00}$ | $121.02_{-17.87}^{+39.66}$ | $59.75_{-8.21}^{+17.93}$ | $-0.03_{-0.19}^{+0.20}$ | $0.12_{-0.21}^{+0.23}$ | $-0.04_{-0.18}^{+0.21}$ | $-2.01_{-0.09}^{+0.31}$ | $0.23_{-0.06}^{+0.22}$
$3.5$ | $-15.14_{-44.56}^{+43.37}$ | $9.16_{-44.03}^{+44.69}$ | $20.22_{-19.01}^{+31.10}$ | $122.44_{-68.48}^{+34.90}$ | $141.26_{-32.38}^{+51.70}$ | $71.17_{-14.71}^{+36.97}$ | $0.04_{-0.63}^{+0.29}$ | $0.11_{-0.33}^{+0.42}$ | $0.01_{-0.25}^{+0.32}$ | $-1.96_{-0.15}^{+0.40}$ | $0.29_{-0.12}^{+0.29}$
$4.0$ | $-17.50_{-47.63}^{+59.88}$ | $14.71_{-54.04}^{+48.99}$ | $22.82_{-24.25}^{+35.65}$ | $126.39_{-56.79}^{+37.60}$ | $153.86_{-41.19}^{+42.42}$ | $77.91_{-17.39}^{+47.86}$ | $0.08_{-0.55}^{+0.30}$ | $0.02_{-0.36}^{+0.33}$ | $0.06_{-0.28}^{+0.31}$ | $-1.81_{-0.28}^{+0.28}$ | $0.40_{-0.20}^{+0.21}$
Without Angular Momentum Constraints
$2.5$ | $-4.59_{-18.19}^{+16.53}$ | $8.02_{-13.30}^{+13.74}$ | $19.82_{-14.06}^{+14.19}$ | $128.47_{-12.97}^{+12.85}$ | $113.62_{-9.52}^{+11.59}$ | $112.02_{-9.69}^{+11.00}$ | $-0.07_{-0.12}^{+0.12}$ | $0.03_{-0.13}^{+0.16}$ | $-0.08_{-0.12}^{+0.12}$ | $-2.04_{-0.06}^{+0.07}$ | $0.21_{-0.04}^{+0.04}$
$3.0$ | $-7.41_{-20.85}^{+19.00}$ | $16.01_{-14.27}^{+15.70}$ | $22.18_{-16.36}^{+16.02}$ | $127.30_{-14.47}^{+14.08}$ | $111.46_{-10.46}^{+13.30}$ | $113.48_{-10.97}^{+12.90}$ | $-0.07_{-0.14}^{+0.14}$ | $0.02_{-0.15}^{+0.16}$ | $-0.10_{-0.13}^{+0.13}$ | $-2.05_{-0.06}^{+0.07}$ | $0.21_{-0.04}^{+0.04}$
$3.5$ | $-12.77_{-25.26}^{+24.33}$ | $18.86_{-19.23}^{+20.68}$ | $19.36_{-18.59}^{+21.11}$ | $132.05_{-17.40}^{+16.77}$ | $112.32_{-12.67}^{+18.44}$ | $117.99_{-13.23}^{+18.75}$ | $-0.01_{-0.17}^{+0.18}$ | $-0.03_{-0.17}^{+0.19}$ | $-0.12_{-0.15}^{+0.15}$ | $-2.08_{-0.07}^{+0.08}$ | $0.20_{-0.04}^{+0.06}$
$4.0$ | $-15.40_{-30.43}^{+28.04}$ | $19.34_{-21.91}^{+25.40}$ | $22.07_{-22.28}^{+24.34}$ | $130.62_{-20.40}^{+20.83}$ | $112.94_{-16.24}^{+19.77}$ | $115.40_{-15.42}^{+18.53}$ | $0.06_{-0.21}^{+0.23}$ | $0.03_{-0.21}^{+0.23}$ | $-0.15_{-0.18}^{+0.17}$ | $-2.10_{-0.07}^{+0.08}$ | $0.19_{-0.04}^{+0.05}$
Figure 15: The corner plot of the model parameters of the disk component for
the sample in the region of $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc. Figure 16:
The corner plot of the model parameters of the substructure component for the
sample in the region of $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc. Figure 17: The
corner plot of the model parameters of the halo component for the sample in
the region of $r\in[7.5,8.5]$ kpc & $|z|>2.5$ kpc.

## References

  * Amsler et al. (2008) Amsler, C., Doser, M., Antonelli, M., et al. 2008, Physics Letters B, 667, 1, doi: https://doi.org/10.1016/j.physletb.2008.07.018
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bahcall & Soneira (1980) Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73, doi: 10.1086/190685
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, The Astrophysical Journal, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, Monthly Notices of the Royal Astronomical Society, 478, 611, doi: 10.1093/mnras/sty982
  * Belokurov & Kravtsov (2022) Belokurov, V., & Kravtsov, A. 2022, Monthly Notices of the Royal Astronomical Society, 514, 689, doi: 10.1093/mnras/stac1267
  * Belokurov et al. (2022) Belokurov, V., Vasiliev, E., Deason, A. J., et al. 2022, Monthly Notices of the Royal Astronomical Society, 518, 6200, doi: 10.1093/mnras/stac3436
  * Bennett & Bovy (2018) Bennett, M., & Bovy, J. 2018, Monthly Notices of the Royal Astronomical Society, 482, 1417, doi: 10.1093/mnras/sty2813
  * Besla et al. (2019) Besla, G., Peter, A., & Garavito-Camargo, N. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 013, doi: 10.1088/1475-7516/2019/11/013
  * Bhattacharjee et al. (2013) Bhattacharjee, P., Chaudhury, S., Kundu, S., & Majumdar, S. 2013, Phys. Rev. D, 87, 083525, doi: 10.1103/PhysRevD.87.083525
  * Bignone et al. (2019) Bignone, L. A., Helmi, A., & Tissera, P. B. 2019, The Astrophysical Journal Letters, 883, L5, doi: 10.3847/2041-8213/ab3e0e
  * Bird et al. (2019) Bird, S. A., Xue, X.-X., Liu, C., et al. 2019, The Astronomical Journal, 157, 104, doi: 10.3847/1538-3881/aafd2e
  * Bird et al. (2021) —. 2021, The Astrophysical Journal, 919, 66, doi: 10.3847/1538-4357/abfa9e
  * Bonaca et al. (2020) Bonaca, A., Conroy, C., Cargile, P. A., et al. 2020, The Astrophysical Journal Letters, 897, L18, doi: 10.3847/2041-8213/ab9caa
  * Bozorgnia et al. (2019) Bozorgnia, N., Fattahi, A., Cerdeño, D. G., et al. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 045, doi: 10.1088/1475-7516/2019/06/045
  * Buder et al. (2021) Buder, S., Lind, K., Ness, M. K., et al. 2021, Monthly Notices of the Royal Astronomical Society, 510, 2407, doi: 10.1093/mnras/stab3504
  * Caldwell & Ostriker (1981) Caldwell, J. A. R., & Ostriker, J. P. 1981, ApJ, 251, 61, doi: 10.1086/159441
  * 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
  * Das et al. (2020) Das, P., Hawkins, K., & Jofré, P. 2020, Monthly Notices of the Royal Astronomical Society, 493, 5195, doi: 10.1093/mnras/stz3537
  * de Salas et al. (2019) de Salas, P., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, Journal of Cosmology and Astroparticle Physics, 2019, 037, doi: 10.1088/1475-7516/2019/10/037
  * Deason et al. (2018) Deason, A. J., Belokurov, V., Koposov, S. E., & Lancaster, L. 2018, The Astrophysical Journal Letters, 862, L1, doi: 10.3847/2041-8213/aad0ee
  * Deason et al. (2019) Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, Monthly Notices of the Royal Astronomical Society, 490, 3426, doi: 10.1093/mnras/stz2793
  * Ding et al. (2021) Ding, P.-J., Xue, X.-X., Yang, C., et al. 2021, AJ, 162, 112, doi: 10.3847/1538-3881/ac0892
  * Donaldson et al. (2022) Donaldson, K., Petersen, M. S., & Peñarrubia, J. 2022, Monthly Notices of the Royal Astronomical Society: Letters, 513, 46, doi: 10.1093/mnrasl/slac031
  * Donlon & Newberg (2023) Donlon, T., & Newberg, H. J. 2023, The Astrophysical Journal, 944, 169, doi: 10.3847/1538-4357/acb150
  * Donlon et al. (2022) Donlon, T., Newberg, H. J., Kim, B., & Lépine, S. 2022, The Astrophysical Journal Letters, 932, L16, doi: 10.3847/2041-8213/ac7531
  * Drimmel & Poggio (2018) Drimmel, R., & Poggio, E. 2018, Research Notes of the AAS, 2, 210, doi: 10.3847/2515-5172/aaef8b
  * Drukier et al. (1986) Drukier, A. K., Freese, K., & Spergel, D. N. 1986, Phys. Rev. D, 33, 3495, doi: 10.1103/PhysRevD.33.3495
  * Eilers et al. (2019) Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, The Astrophysical Journal, 871, 120, doi: 10.3847/1538-4357/aaf648
  * Evans et al. (2019) Evans, N. W., O’Hare, C. A. J., & McCabe, C. 2019, Phys. Rev. D, 99, 023012, doi: 10.1103/PhysRevD.99.023012
  * Evans et al. (2016) Evans, N. W., Sanders, J. L., Williams, A. A., et al. 2016, Monthly Notices of the Royal Astronomical Society, 456, 4506, doi: 10.1093/mnras/stv2729
  * Fairbairn et al. (2013) Fairbairn, M., Douce, T., & Swift, J. 2013, Astroparticle Physics, 47, 45, doi: https://doi.org/10.1016/j.astropartphys.2013.06.003
  * Fattahi et al. (2019) Fattahi, A., Belokurov, V., Deason, A. J., et al. 2019, Monthly Notices of the Royal Astronomical Society, 484, 4471, doi: 10.1093/mnras/stz159
  * Feuillet et al. (2020) Feuillet, D. K., Feltzing, S., Sahlholdt, C. L., & Casagrande, L. 2020, Monthly Notices of the Royal Astronomical Society, 497, 109, doi: 10.1093/mnras/staa1888
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Frandsen et al. (2012) Frandsen, M. T., Kahlhoefer, F., McCabe, C., Sarkar, S., & Schmidt-Hoberg, K. 2012, Journal of Cosmology and Astroparticle Physics, 2012, 024, doi: 10.1088/1475-7516/2012/01/024
  * Freese et al. (1988) Freese, K., Frieman, J., & Gould, A. 1988, Phys. Rev. D, 37, 3388, doi: 10.1103/PhysRevD.37.3388
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211.  https://arxiv.org/abs/2208.00211
  * Gallart et al. (2019) Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, Nature Astronomy, 3, 932, doi: 10.1038/s41550-019-0829-5
  * Garrison-Kimmel et al. (2016) Garrison-Kimmel, S., Bullock, J. S., Boylan-Kolchin, M., & Bardwell, E. 2016, Monthly Notices of the Royal Astronomical Society, 464, 3108, doi: 10.1093/mnras/stw2564
  * Gómez et al. (2010) Gómez, F. A., Helmi, A., Brown, A. G. A., & Li, Y.-S. 2010, Monthly Notices of the Royal Astronomical Society, 408, 935, doi: 10.1111/j.1365-2966.2010.17225.x
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 618, L10, doi: 10.1051/0004-6361/201834294
  * Green (2017) Green, A. M. 2017, Journal of Physics G: Nuclear and Particle Physics, 44, 084001, doi: 10.1088/1361-6471/aa7819
  * Guo et al. (2020) Guo, R., Liu, C., Mao, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 495, 4828, doi: 10.1093/mnras/staa1483
  * Guo et al. (2022) Guo, R., Shen, J., Li, Z.-Y., Liu, C., & Mao, S. 2022, The Astrophysical Journal, 936, 103, doi: 10.3847/1538-4357/ac86cd
  * Hagen, Jorrit H. J. & Helmi, Amina (2018) Hagen, Jorrit H. J., & Helmi, Amina. 2018, A&A, 615, A99, doi: 10.1051/0004-6361/201832903
  * Hattori et al. (2021) Hattori, K., Valluri, M., & Vasiliev, E. 2021, Monthly Notices of the Royal Astronomical Society, 508, 5468, doi: 10.1093/mnras/stab2898
  * Helm (1956) Helm, R. H. 1956, Phys. Rev., 104, 1466, doi: 10.1103/PhysRev.104.1466
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85, doi: 10.1038/s41586-018-0625-x
  * Helmi & White (1999) Helmi, A., & White, S. D. M. 1999, Monthly Notices of the Royal Astronomical Society, 307, 495, doi: 10.1046/j.1365-8711.1999.02616.x
  * Herzog-Arbeitman et al. (2018a) Herzog-Arbeitman, J., Lisanti, M., Madau, P., & Necib, L. 2018a, Phys. Rev. Lett., 120, 041102, doi: 10.1103/PhysRevLett.120.041102
  * Herzog-Arbeitman et al. (2018b) Herzog-Arbeitman, J., Lisanti, M., & Necib, L. 2018b, Journal of Cosmology and Astroparticle Physics, 2018, 052, doi: 10.1088/1475-7516/2018/04/052
  * Huang et al. (2016) Huang, Y., Liu, X.-W., Yuan, H.-B., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 2623, doi: 10.1093/mnras/stw2096
  * Ibata et al. (2021) Ibata, R., Malhan, K., Martin, N., et al. 2021, The Astrophysical Journal, 914, 123, doi: 10.3847/1538-4357/abfcc2
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, The Astrophysical Journal, 779, 102, doi: 10.1088/0004-637X/779/2/102
  * Koppelman et al. (2020) Koppelman, H. H., Bos, R. O. Y., & Helmi, A. 2020, A&A, 642, L18, doi: 10.1051/0004-6361/202038652
  * Lancaster et al. (2019) Lancaster, L., Koposov, S. E., Belokurov, V., Evans, N. W., & Deason, A. J. 2019, Monthly Notices of the Royal Astronomical Society, 486, 378, doi: 10.1093/mnras/stz853
  * Lewin & Smith (1996) Lewin, J., & Smith, P. 1996, Astroparticle Physics, 6, 87, doi: https://doi.org/10.1016/S0927-6505(96)00047-3
  * Lisanti et al. (2015) Lisanti, M., Spergel, D. N., & Madau, P. 2015, The Astrophysical Journal, 807, 14, doi: 10.1088/0004-637X/807/1/14
  * Liu et al. (2014) Liu, C., Deng, L.-C., Carlin, J. L., et al. 2014, The Astrophysical Journal, 790, 110, doi: 10.1088/0004-637x/790/2/110
  * 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. (2018) Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2018, Monthly Notices of the Royal Astronomical Society, 482, 3426, doi: 10.1093/mnras/sty2955
  * Malhan & Ibata (2019) Malhan, K., & Ibata, R. A. 2019, Monthly Notices of the Royal Astronomical Society, 486, 2995, doi: 10.1093/mnras/stz1035
  * Meng et al. (2021) Meng, Y., Wang, Z., Tao, Y., et al. 2021, Phys. Rev. Lett., 127, 261802, doi: 10.1103/PhysRevLett.127.261802
  * Myeong et al. (2018) Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018, The Astrophysical Journal Letters, 863, L28, doi: 10.3847/2041-8213/aad7f7
  * Naidu et al. (2020) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, The Astrophysical Journal, 901, 48, doi: 10.3847/1538-4357/abaef4
  * Necib et al. (2019a) Necib, L., Lisanti, M., & Belokurov, V. 2019a, The Astrophysical Journal, 874, 3, doi: 10.3847/1538-4357/ab095b
  * Necib et al. (2019b) Necib, L., Lisanti, M., Garrison-Kimmel, S., et al. 2019b, The Astrophysical Journal, 883, 27, doi: 10.3847/1538-4357/ab3afc
  * Necib et al. (2020) Necib, L., Ostdiek, B., Lisanti, M., et al. 2020, The Astrophysical Journal, 903, 25, doi: 10.3847/1538-4357/abb814
  * Nelson et al. (2019) Nelson, D., Pillepich, A., Springel, V., et al. 2019, Monthly Notices of the Royal Astronomical Society, 490, 3234, doi: 10.1093/mnras/stz2306
  * O’Hare et al. (2020) O’Hare, C. A. J., Evans, N. W., McCabe, C., Myeong, G., & Belokurov, V. 2020, Phys. Rev. D, 101, 023006, doi: 10.1103/PhysRevD.101.023006
  * Ostriker et al. (1974) Ostriker, J. P., Peebles, P. J. E., & Yahil, A. 1974, ApJ, 193, L1, doi: 10.1086/181617
  * Pila-Díez, B. et al. (2015) Pila-Díez, B., de Jong, J. T. A., Kuijken, K., van der Burg, R. F. J., & Hoekstra, H. 2015, A&A, 579, A38, doi: 10.1051/0004-6361/201425457
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, Monthly Notices of the Royal Astronomical Society, 490, 3196, doi: 10.1093/mnras/stz2338
  * Price-Whelan (2018) Price-Whelan, A. 2018, adrn/pyia: v0.2, v0.2, Zenodo, doi: 10.5281/zenodo.1228136
  * Reid & Brunthaler (2004) Reid, M. J., & Brunthaler, A. 2004, The Astrophysical Journal, 616, 872, doi: 10.1086/424960
  * Sahlholdt et al. (2019) Sahlholdt, C. L., Casagrande, L., & Feltzing, S. 2019, The Astrophysical Journal Letters, 881, L10, doi: 10.3847/2041-8213/ab321e
  * Schönrich (2012) Schönrich, R. 2012, Monthly Notices of the Royal Astronomical Society, 427, 274, doi: 10.1111/j.1365-2966.2012.21631.x
  * Schönrich & Binney (2011) Schönrich, R., & Binney, J. 2011, Monthly Notices of the Royal Astronomical Society, 419, 1546, doi: 10.1111/j.1365-2966.2011.19816.x
  * Schwarz (1978) Schwarz, G. 1978, The Annals of Statistics, 6, 461 , doi: 10.1214/aos/1176344136
  * Sesar et al. (2011) Sesar, B., Jurić, M., & Ivezić, Ž. 2011, The Astrophysical Journal, 731, 4, doi: 10.1088/0004-637X/731/1/4
  * Smith et al. (2009) Smith, M. C., Evans, N. W., & An, J. H. 2009, The Astrophysical Journal, 698, 1110, doi: 10.1088/0004-637X/698/2/1110
  * Smith-Orlik et al. (2023) Smith-Orlik, A., Ronaghi, N., Bozorgnia, N., et al. 2023, The impact of the Large Magellanic Cloud on dark matter direct detection signals.  https://arxiv.org/abs/2302.04281
  * Tian et al. (2015) Tian, H.-J., Liu, C., Carlin, J. L., et al. 2015, The Astrophysical Journal, 809, 145, doi: 10.1088/0004-637X/809/2/145
  * Vasiliev et al. (2020) Vasiliev, E., Belokurov, V., & Erkal, D. 2020, Monthly Notices of the Royal Astronomical Society, 501, 2279, doi: 10.1093/mnras/staa3673
  * Vincenzo et al. (2019) Vincenzo, F., Spitoni, E., Calura, F., et al. 2019, Monthly Notices of the Royal Astronomical Society: Letters, 487, L47, doi: 10.1093/mnrasl/slz070
  * Wang et al. (2011) Wang, J., Navarro, J. F., Frenk, C. S., et al. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1373, doi: 10.1111/j.1365-2966.2011.18220.x
  * Watkins et al. (2009) Watkins, L. L., Evans, N. W., Belokurov, V., et al. 2009, Monthly Notices of the Royal Astronomical Society, 398, 1757, doi: 10.1111/j.1365-2966.2009.15242.x
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, Monthly Notices of the Royal Astronomical Society, 183, 341, doi: 10.1093/mnras/183.3.341
  * Wu et al. (2022) Wu, W., Zhao, G., Xue, X.-X., Bird, S. A., & Yang, C. 2022, ApJ, 924, 23, doi: 10.3847/1538-4357/ac31ac
  * Wu et al. (2011) Wu, Y., Luo, A. L., Li, H.-N., et al. 2011, Research in Astronomy and Astrophysics, 11, 924, doi: 10.1088/1674-4527/11/8/006
  * Xia et al. (2016) Xia, Q., Liu, C., Mao, S., et al. 2016, Monthly Notices of the Royal Astronomical Society, 458, 3839, doi: 10.1093/mnras/stw565
  * Xue et al. (2015) Xue, X.-X., Rix, H.-W., Ma, Z., et al. 2015, The Astrophysical Journal, 809, 144, doi: 10.1088/0004-637X/809/2/144
  * Xue et al. (2014) Xue, X.-X., Ma, Z., Rix, H.-W., et al. 2014, The Astrophysical Journal, 784, 170, doi: 10.1088/0004-637x/784/2/170
  * Yang et al. (2019) Yang, C., Xue, X.-X., Li, J., et al. 2019, The Astrophysical Journal, 886, 154, doi: 10.3847/1538-4357/ab48e2
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723, doi: 10.1088/1674-4527/12/7/002