# The Kinematic Properties of Milky Way Stellar Halo Populations

James M. M. Lane1, Jo Bovy1, & J. Ted Mackereth1,2,3
1Department of Astronomy and Astrophysics, University of Toronto, 50 St.
George Street, Toronto ON, M5S 3H4, Canada
2Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto ON, M5S 3H8, Canada
3Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto ON, M5S 3H4, Canada E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In the Gaia era stellar kinematics are extensively used to study Galactic halo
stellar populations, to search for halo structures, and to characterize the
interface between the halo and hot disc populations. We use distribution
function-based models of modern datasets with 6D phase space data to
qualitatively describe a variety of kinematic spaces commonly used in the
study of the Galactic halo. Furthermore, we quantitatively assess how well
each kinematic space can separate radially anisotropic from isotropic halo
populations. We find that scaled action space (the “action diamond”) is
superior to other commonly used kinematic spaces at this task. We present a
new, easy to implement selection criterion for members of the radially-
anisotropic Gaia-Enceladus merger remnant, which we find achieves a sample
purity of 82 per cent in our models with respect to contamination from the
more isotropic halo. We compare this criterion to literature criteria, finding
that it produces the highest purity in the resulting samples, at the expense
of a modest reduction in completeness. We also show that selection biases that
underlie nearly all contemporary spectroscopic datasets can noticeably impact
the $E-L_{z}$ distribution of samples in a manner that may be confused for
real substructure. We conclude by providing recommendations for how authors
should use stellar kinematics in the future to study the Galactic stellar
halo.

###### keywords:

Galaxy: halo – Galaxy: kinematics and dynamics – Galaxy: solar neighbourhood –
Galaxy: stellar content

††pubyear: 2021††pagerange: The Kinematic Properties of Milky Way Stellar Halo
Populations–The Kinematic Properties of Milky Way Stellar Halo Populations

##  1 Introduction

It is an auspicious time to study the dynamics and structure of our Milky Way
galaxy, in no small part thanks to a recent wealth of data obtained through
large astronomical surveys. Chief among these is the astrometric space mission
Gaia (Gaia Collaboration et al., 2016) which, as part of its most recent data
release, has provided 5-d kinematic phase space information for nearly 1.5
billion stars, and full 6-d phase space information for a further $\sim 7$
million stars (Gaia Collaboration et al., 2020). In parallel to this project,
massive spectroscopic campaigns such as APOGEE (Majewski et al., 2017), LAMOST
(Cui et al., 2012), Gaia-ESO (Gilmore et al., 2012), GALAH (Martell et al.,
2017), SEGUE (Yanny et al., 2009), RAVE (Steinmetz et al., 2020, for the most
recent data release), and H3 (Conroy et al., 2019) have provided accurate
radial velocities, chemical abundances, and age estimates for hundreds of
thousands of these stars. Combined, these two efforts have yielded a rich,
multidimensional data set which has not only enabled novel discoveries, but in
concert with powerful statistical techniques will allow us to discern among
competing models of the structure, history, and composition of each part of
our Milky Way galaxy with exquisite detail.

The nature of the stellar halo, in particular, has come into much sharper
focus in the Gaia era. Since Gaia’s second data release, a significant amount
of evidence has been amassed that the Milky Way was subject to a handful of
distinct, significant mergers early in its life. The most well-studied example
to date is the colourfully named “Sausage”, or “Gaia-Enceladus” (GE), which is
a population of metal-rich halo stars with uncharacteristically large radial
motions (Helmi et al., 2018; Belokurov et al., 2018; Haywood et al., 2018).
This overdensity stands out in phase space from the remainder of halo stars,
which exhibit more isotropically-distributed motions. Estimates of the mass of
this population are in the range of $10^{8}$ to $10^{10}$ M⊙(Belokurov et al.,
2018; Helmi et al., 2018; Fattahi et al., 2019; Deason et al., 2019; Mackereth
et al., 2019a; Mackereth & Bovy, 2020; Kruijssen et al., 2020; Naidu et al.,
2021) which would place its stellar mass among the contemporary big satellites
of the Milky Way, such as Sagittarius and the Large Magellanic Cloud (e.g.
Niederste-Ostholt et al., 2010; van der Marel, 2006).

GE is only one among an assortment of spatially incoherent halo structures
which have either been discovered or extensively studied thanks to Gaia
kinematics and supporting spectroscopy. These include Thamnos (Koppelman et
al., 2019b), the Sequoia (Myeong et al., 2019), the Helmi streams (Helmi et
al., 1999; Koppelman et al., 2018, 2019a), the Inner Galaxy Structure (IGS,
also known as Heracles; Horta et al., 2021), as well as numerous structures
identified by Naidu et al. (2020). Paralleling these results driven by stellar
data, examination of the Galactic globular clusters population points to a
merger known as “Kraken” (Kruijssen et al., 2020) or “Koala” (Forbes, 2020)
which predates the above-mentioned events. While many of these structures
require further study to confirm their nature and origin, it is nontheless
clear that significant merger events play a major role in shaping the stellar
halo we see today.

The stellar halo is not solely composed of ex-situ accretion remnants, but
also in-situ stars formed within the Milky Way. Mounting kinematic and
chemical evidence suggests that this in-situ component consists mostly of
stars from the thicker components of the disc (Helmi et al., 2018; Haywood et
al., 2018; Di Matteo et al., 2019). This idea, while not new (e.g. Nissen &
Schuster, 2010), is put into new perspective given the emerging importance of
mergers in the early evolution of the Milky Way. Additionally, an emerging
class of halo substructures known as “Splash” populations, which may be
remnants of a proto-disc disrupted by an infalling satellite, may be evidence
of the dynamical link between halo and thick disc (Belokurov et al., 2020;
Naidu et al., 2020; Lian et al., 2020).

While each halo population is often discussed as being “distinct” in terms of
its spatial, kinematic, and/or chemical properties, there is significant
overlap in the properties of many, if not all, of these structures. This
presents an issue when attempting to dissect the stellar halo on a granular
level, as is the intent with current and upcoming studies. We recognize a gap
in the established literature, whereby there does not exist simple kinematic
models to contextualize these large range of data. N-body models are commonly
used to interpret results, but tailoring the simulations to individual
datasets is challenging.

In this work we examine mock samples of stars generated from distribution
functions which represent the disc and halo populations commonly studied in
the Milky Way. Our goals are twofold: first we seek to provide context for
observed kinematics from the current and next generation of surveys. Second we
use these models to numerically evaluate the best way to separate radially
biased halo populations from more isotropic halo populations, an effective
proxy for isolating GE remnant stars.

##  2 Model Potential and Distribution Functions

Here we introduce the distribution functions (DFs) and accompanying
gravitational potentials that we use to model the Milky Way disc and halo.
First, we present the potential and then we describe the DFs and the
techniques used to draw kinematic samples from them. Throughout we assume
that: the Sun is located at $(R,z,\phi)$ =
$(8.178~{}\text{kpc},0.028~{}\text{kpc},0)$ (Gravity Collaboration et al.,
2019; Bennett & Bovy, 2019), and the solar motion is
$(U,V,W)=(11.1,12.24,7.25)$ km s-1 (Schönrich et al., 2010). We use a left-
handed Galactocentric cylindrical coordinate system where angular momentum is
positive in the direction of Galactic rotation. By default we express
quantities in cylindrical coordinates and are explicit when we do otherwise.

###  2.1 Milky Way Potential

To represent the Milky Way we use the MWPotential2014 model from Bovy (2015).
The potential has three components. The disc is represented by a Miyamoto &
Nagai (1975) profile with mass $M_{d}=6.8\times 10^{10}$ M⊙, a scale length of
$3$ kpc, and a scale height of $280$ pc. The bulge is an exponentially-
truncated power law density profile with mass $M_{b}=0.5\times 10^{10}$ M⊙,
power law index $-1.8$, and exponential scale radius of $1.9$ kpc. The halo is
an NFW (Navarro et al., 1997) profile with a virial mass of $8\times 10^{11}$
M⊙ and concentration of 15.3. For a discussion of the derivation of the model
parameters which describe MWPotential2014, see section 3.5 of Bovy (2015).

###  2.2 Disc Distribution Functions

We approximate the disc’s intricate structure (e.g., Bovy et al. 2012) as a
simple “thin” and “thick” superposition and to represent both disc components
we use the quasi-isothermal DF of Binney & McMillan (2011) (see also Binney,
2010). This DF takes the form

$f(J_{r},L_{z},J_{z})=f_{\sigma_{r}}(J_{r},L_{z})f_{\sigma_{z}}(J_{z})\,,$ (1)

where $f_{\sigma_{r}}$ largely controls the dynamics within the plane of the
disc, and $f_{\sigma_{z}}$ largely controls the vertical dynamics
perpendicular to the disc. Each of these are functions of the actions $J_{r}$,
$L_{z}$, and $J_{z}$. $f_{\sigma_{r}}$ is given by

$f_{\sigma_{r}}=\frac{\Omega\Sigma}{\pi\sigma_{r}^{2}\kappa}\bigg{\rvert}_{\small{R_{c}}}[1+\tanh(L_{z}/L_{0})]\mathrm{e}^{-\kappa
J_{r}/\sigma_{r}^{2}}\,,$ (2)

where $\Sigma$ and $\Omega$ are the local surface density and orbital
frequency respectively. $R_{c}$ and $\kappa$ are the guiding center radius and
epicycle frequency of a given orbit. $L_{0}$ is a cutoff angular momentum.
$\sigma_{r}$ is the radial velocity dispersion profile. $\Sigma$ has the
exponential form

$\Sigma(L_{z})=\Sigma_{0}\mathrm{e}^{(R_{0}-R_{c})/R_{d}}$ (3)

where $\Sigma_{0}$ is an arbitrary normalization, $R_{0}$ is a scale radius
(Which we set to the location of the Sun), and $R_{d}$ is the scale length of
the disc. The second component, $f_{\sigma_{z}}$, is given by a pseudo-
isothermal DF

$f_{\sigma_{z}}=\frac{\nu_{z}}{2\pi\sigma_{z}^{2}}\mathrm{e}^{-\nu_{z}J_{z}/\sigma_{z}^{2}}\,,$
(4)

where $\nu_{z}$ and $J_{z}$ are the vertical epicycle frequency and action and
$\sigma_{z}$ is the vertical velocity dispersion profile. The parameters we
use for the thin and thick disc DFs are summarized in Table 1. Both the radial
and vertical exponential profiles have exponential forms given by111Binney &
McMillan (2011) have a multiplicative factor $q$ in these exponentials, but it
is unity here.

$\begin{split}\sigma_{r}(L_{z})=&\sigma_{r,0}\mathrm{e}^{(R_{0}-R_{c})/R_{d}}\\\
\sigma_{r}(L_{z})=&\sigma_{z,0}\mathrm{e}^{(R_{0}-R_{c})/R_{d}}\,.\\\
\end{split}$

Here $\sigma_{[r,z],0}$ are characteristic velocity dispersions at $R_{0}$. We
use, overall, hotter disc DF parameters than are advocated for by Binney &
McMillan (2011), but that are not unlike other measurements have suggested
(e.g. Bovy et al., 2012; Mackereth et al., 2019b). By choosing disc DF
parameters on the hot end of what is normally used we hedge our results
towards including disc samples which intrude most on the kinematic space
normally assigned to the halo.

Table 1: Properties of quasi-isothermal thin and thick disc DFs Parameter | Thin disc | Thick disc
---|---|---
$R_{d}$ | 3 kpc | 2 kpc
$R_{0}$ | 8.178 kpc | 8.178 kpc
$\sigma_{r,0}$ | 40 km s-1 | 60 km s-1
$\sigma_{z,0}$ | 30 km s-1 | 45 km s-1
$L_{0}$ | 10 kpc km s-1 | 10 kpc km s-1

To sample from this DF, we use galpy222https://github.com/jobovy/galpy (Bovy,
2015) where it is implemented as the class galpy.df.quasiisothermaldf. The
phase space samples we generate are described further on in § 3, but for our
purposes we are only interested in sampling velocities from this DF given
fixed positions. Velocity samples are obtained using rejection sampling of the
DF at each individual radius and height above the midplane. Actions for the DF
are calculated using the “Stäckel fudge” method of Binney (2012). In this
approach our Milky Way potential is locally approximated for each individual
sample orbit as a Stäckel potential using a focal length estimated with
equation (9) of Sanders (2012).

###  2.3 The Stellar Halo Distribution Function

For the stellar halo we use spherical DFs with constant anisotropy. These
types of DFs generally take the form

$f(\mathcal{E},L)=L^{-2\beta}f_{1}(\mathcal{E})\,.$ (5)

Conventionally, here, $\mathcal{E}=\Psi-\frac{1}{2}v^{2}$ is the binding
energy, where $\Psi=-\Phi+\Phi_{0}$ is the negative of the gravitational
potential offset to be 0 at infinity, and $v$ is the magnitude of the
velocity. $L$ is the total 3-dimensional angular momentum. The parameter
$\beta$ is the orbital anisotropy, defined as

$\beta=1-\frac{\sigma^{2}_{\theta}+\sigma^{2}_{\phi}}{2\sigma^{2}_{r}}\,,$ (6)

where $\sigma_{[\theta,\phi,r]}$ are the spherical velocity dispersions in the
polar, azimuthal, and radial directions respectively. $\beta$ ranges from
$-\infty$ (completely tangential orbits), over 0 (ergodic, isotropic, or
unbiased orbits), and up to 1 (completely radial orbits). As shown in chapter
4.3 of Binney & Tremaine (2008) the function $f_{1}$ is related to the mass
density of the DF333Note they use number density instead of mass density, but
that just amounts to an overall normalization of the DF, $\rho$ as

$\frac{2^{\beta-1/2}}{2\pi
I_{\beta}}r^{2\beta}\rho=\int_{0}^{\Psi}\mathrm{d}\mathcal{E}\frac{f_{1}(\mathcal{E})}{(\Psi-\mathcal{E})^{\beta-1/2}}\,.$
(7)

$I_{\beta}$ is a constant given by

$I_{\beta}=\sqrt{\pi}\frac{\Gamma(1-\beta)}{\Gamma(3/2-\beta)}\,,$ (8)

where $\Gamma$ is the usual Gamma function.

For $1/2<\beta<1$, a DF of this form is an Abel integral equation that can be
inverted to yield, often numerically, a solution for $f_{1}(\mathcal{E})$. For
lower values of $\beta$, one or more derivatives with respect to $\Psi$ turn
this expression into an invertible Abel integral equation (Cuddeford, 1991).
Note that while we can use this DF for $\beta>1/2$, generally it has an
increasingly large region where the DF is negative and, thus, unphysical. We
set the DF to zero wherever the expression gives a negative number. This
causes a slight, but unimportant inconsistency between the DF and the density
profile. For $\beta=0.9$ (the highest value we will consider) the profile
becomes shallower, with density decreasing by about 10 per cent at 3 kpc being
roughly unchanged at the solar circle, and increasing by a few per cent at
more distant radii.

For many implementations of this DF where the system is self-gravitating the
potential, $\Psi$, and the density profile, $\rho$, are a potential-density
pair. This is not necessary, however, and $\rho$ may be independent from
$\Psi$. For our DFs $\rho$ is an exponentially truncated power law density
profile of the form

$\rho(r)\propto r^{-\alpha}\mathrm{e}^{r/r_{c}},$ (9)

where $\alpha$ is the inner power law slope and $r_{c}$ is the truncation
radius. The overall normalization of the density profile is unimportant
because it only determines positions and kinematics of samples in a relative
sense. We set $\alpha=3.5$ and $r_{c}=30$ kpc. Our choice of parameters is
inspired by Mackereth & Bovy (2020), who find $\alpha=3.49$ and $r_{c}=25$ kpc
fit a large sample of halo stars spanning both [Fe/H] and [$\alpha$/Fe]. These
fits stand in broad agreement with prior results, both in terms of the power
law slope, and the transition radius at which the density profile steepens
(Deason et al., 2011; Xue et al., 2015; Iorio et al., 2018).

The potential in the calculation of these stellar-halo DFs is spherical.
Because our adopted MWPotential2014 model is non-spherical, we construct a
spherically symmetric analog of this potential by calculating its enclosed
mass within spherical shells as a function of Galactocentric radius $M(<r)$.
The spherical version of the potential is then defined as
$\Phi_{MW}=-\frac{GM(<r)}{r}$. While possible, it is difficult to create
anisotropic DFs for spherical systems, such as the stellar halo, with a given
density and anisotropy embedded within axisymmetric, non-spherical potentials
(Posti et al., 2015; Williams & Evans, 2015). Our approach, while approximate,
is simple to implement and gives us greater control over the kinematic and
spatial properties of our halo DFs. The act of sphericalizing the potential
should have a minimal impact on our results, which focuses on kinematic
properties which are sensitive to the total enclosed mass of the potential
rather than increased mass near the disc plane.

We sample from DFs of this type using galpy, where they are implemented as the
class galpy.df.constantbetadf. More details on our implementation of these
spherical DFs and how we specifically sample from them can be found in Bovy &
Lane (2021; in preparation).

##  3 DF Samples

In this section, we present the samples used to explore the observable
kinematics of Milky Way stellar populations, and describe how we calculate
their kinematics. We first outline the overall properties of our samples and
then describe how we generate the two specific sets of samples that we use
throughout the rest of the paper. Finally, we present the kinematics of each
of these samples.

Broadly speaking our samples consist of phase space positions drawn from a
combination of the distribution functions outlined in § 2.2 and 2.3,
representing the disc (approximated as a thin plus thick disc) and stellar
halo. We do not consider the stellar bulge for two reasons. First, the bulge
does not contribute stellar density near the solar neighbourhood in any
meaningful capacity, and therefore bulge stars are reasonably easy to remove
from real stellar samples by cutting stars in angular proximity to the
Galactic centre. This means one can be confident that the bulge does not
represent a substantial source of contamination when studying halo
populations, in contrast to the stellar disc. Second, the complexity of bulge
kinematics due, for example, to the presence of the rotating Galactic bar
means that it is challenging to realistically model bulge populations using
DFs.

###  3.1 The _Solar_ Samples

The first set of samples, which we refer to as _Solar_ are located at the
Sun’s position: ($R,z$) = (8.178 kpc, 20.8 pc) (Gravity Collaboration et al.,
2019; Bennett & Bovy, 2019). This set is designed to be instructive, and so
having all samples share the same position simplifies the kinematics. Yet,
because spatial gradients in the DF occur over multiple-kpc scales, this
sample is also somewhat realistic in that it represents stellar samples
confined to the $\approx 1$ kpc region around the Sun, a distance over which
the halo DF varies little. This includes existing samples such as the Gaia
radial velocity sample (see section 3.3 in Gaia Collaboration et al., 2020,
for the most recent data) or GALAH (Martell et al., 2017; Buder et al., 2020).
The _Solar_ set of samples consists of six independent samples of stars, which
differ by the value of $\beta$ used for the stellar halo DF.

Each of the _Solar_ samples consists of $10^{4}$ thin disc orbits and
$1.2\times 10^{3}$ thick disc orbits, sampled according to the procedure
outlined in § 2.2. These numbers are chosen such that the ratio of thick to
thin disc stars is approximately 1:8.5, which reflects the ratio of thick to
thin disc stars expected at the solar circle (e.g. see the potential of
McMillan, 2017). Each of the samples also contains $10^{4}$ orbits sampled
from our spherical halo DF with a different value of $\beta$. We are not
concerned with the relative numbers of halo and disc stars, because we will
differentiate between the two populations when studying their properties. The
first five samples each use a different value of $\beta$ from among
$[0,-0.5,0.5,0.7,0.9]$. The sixth sample has a composite stellar halo, also
$10^{4}$ orbits total, consisting of a 50:50 mixture of orbits from DFs with
$\beta=0.5$ and 0.9 respectively. We refer to this final sample as our
fiducial model.

The ratio of the component DFs in this model follows from the recent Milky Way
stellar halo mass measurements of Mackereth & Bovy (2020), who find that the
ratio of mass in low and high eccentricity stars is approximately 50:50. These
low and high eccentricity populations are interpreted as a traditional metal
poor stellar halo with $\beta\approx 0.2-0.4$, and the more metal rich Gaia-
Enceladus remnant which has been shown to be characterized by large
anisotropies of $\beta$ close to 0.9 (Belokurov et al., 2018; Fattahi et al.,
2019). This model is broadly consistent with other results that have suggested
that these two halo components exist in roughly equal proportion near the
solar neighbourhood (Belokurov et al., 2018; Lancaster et al., 2019). By using
$\beta=0.5$ to represent the more isotropic component of the halo we err on
the side of slightly larger $\beta$, but this allows us to be confident we do
not underestimate the degree of contamination from this component when
analyzing their kinematics. This is similar to the reasoning behind our choice
of slightly hotter disc DF parameters discussed in § 2.2.

Throughout the paper we will refer to the two constituent components of our
fiducial model as “low-$\beta$” and “high-$\beta$” for the $\beta=0.5$ and
$\beta=0.9$ halos respectively. While our fiducial model is inspired by
literature results which employ metallicity and eccentricity as discriminators
for the halo populations, we emphasize that we define our halo components
using $\beta$, and therefore they may not be perfectly faithful
representations of a halo defined in another manner.

###  3.2 The _Survey_ Samples

The second set of samples has physical positions which are matched to the
locations of stars in the APOGEE DR16 catalogue (Ahumada et al., 2020), rather
than occupying a single location. We refer to these as the _Survey_ samples.
These samples more realistically represent modern and upcoming datasets where
Gaia data is augmented by spectroscopy, which allows for the calculation of
accurate spectro-photometric parallaxes (e.g., Leung & Bovy, 2019b) and radial
velocities for a small number of stars (compared to the broader Gaia dataset)
out to large ($\sim 20$ kpc) distances. Examples of these surveys include the
contemporary APOGEE Ahumada et al. (2020) and H3 (Conroy et al., 2019)
surveys, as well as the future WEAVE (Dalton et al., 2016), 4MOST (de Jong et
al., 2019), DESI (Levi et al., 2019) and SDSS-V (Kollmeier et al., 2017)
projects.

To generate the _Survey_ sample we start with the on-sky locations, [Fe/H] and
[Mg/Fe] abundances of stars in the APOGEE DR16 release (Ahumada et al., 2020).
We supplement this data with accurate spectro-photometric distances derived
using the Bayesian neural network framework astroNN (Leung & Bovy, 2019a, b)
(note these are based on Gaia DR2 (Gaia Collaboration et al., 2018)). Our goal
is to use this information to obtain a rough set of positions corresponding to
three populations we seek to study: the halo, thin disc, and thick disc.

We perform a number of cuts on the APOGEE data to ensure both quality and that
the underlying stellar populations correctly reflect our model DFs. The cuts
are as follows. The $\log g$ uncertainty is less than 0.1, which removes many
dwarfs with poorly constrained properties. The relative distance uncertainty,
$\sigma_{d}/d$, is less than 0.2. We remove fields near the Galactic bulge,
which we define as those with $-20^{\circ}<\ell<20^{\circ}$ and
$|b|<20^{\circ}$. We also remove any fields which contain distinct groups of
stars from a globular cluster. We first identify which APOGEE fields contain
the on-sky location of a known globular cluster from the catalogue of Harris
(1996, The December 2010 version). We then examine the distribution of [Fe/H]
abundances and radial velocities for the stars in the field, in particular
those within two tidal radii of the globular cluster location. Fields where
there is an obvious enhancement of stars with the appropriate [Fe/H] and
radial velocities near the on-sky location of the globular cluster are removed
from consideration. Figure 1 shows the spatial distribution of the final
selection of APOGEE stars. Our sample spans a wide range of the Milky Way,
nearly 15 kpc in each direction from the location of the Sun.

Figure 1: The APOGEE DR16 sample shown in Galactocentric coordinates. The top
panel shows Galactocentric cylindrical radius and height above the disc, the
bottom panel shows Galactocentric X and Y, where X is along the Sun-Galactic
Centre line, and Y is positive in the direction of Galactic rotation. In both
panels the Sun is shown at (X,Y,Z) = (8.178 kpc, 0, 0.028 kpc).

We separate the APOGEE stars into three categories: halo, thin disc, and thick
disc using their [Fe/H] and [Mg/Fe] abundances as well as their velocities. We
stress that the division into a thin and thick disc is purely an approximation
to the more complex structure of the disc, but one that is good enough for our
main purpose of investigating the stellar halo. Figure 2 shows these
abundances for our final APOGEE sample, as well as the boundaries we adopt to
separate each stellar population. We choose these boundaries based on the
number density of APOGEE stars in the abundance space, as well as their
eccentricities (which color-codes the points in the top panel of Figure 2).
The thin-thick disc boundary follows the upper edge of the low-[Mg/Fe] disc
locus, while the disc-halo boundary follows the low-[Fe/H] edge of the disc
locus and approximately matches the transition from low to high eccentricity.
We also add to the halo sample any stars lying in the thin or thick disc
chemical regions which have net velocities with respect to their local
circular velocity frame which are greater than the magnitude of their local
circular velocity:
$\lvert\vec{v}-\vec{v}_{\mathrm{circ}}(R)\rvert>\lvert\vec{v}_{\mathrm{circ}}(R)\rvert$.
We find that this efficiently attributes many higher metallicity stars with
genuine halo kinematics to the halo sample. We examine the energy and angular
momentum plane, as well as the tangential and radial velocity plane for each
of these selections to ensure that they reflect what we expect from the
kinematics of the respective populations. Overall, however, the accuracy of
this separation is not a major concern, because this only sets the positions
based on which velocities will be sampled from the individual DFs.

The final separated samples consist of 155,935 thin disc stars, 29,634 thick
disc stars, and 2,798 halo stars. Like the _Solar_ sample, the _Survey_
samples consist of six individual sets of orbits which differ by the value of
$\beta$ used for the halo DF. Again, the first five samples each use a
different value of $\beta$ from among $[0,-0.5,0.5,0.7,0.9]$, while the sixth
sample has a composite stellar halo consisting of a 50:50 mixture of orbits
from DFs with $\beta=0.5$ and $0.9$ respectively.

Figure 2: [Mg/Fe] versus [Fe/H] abundances for the APOGEE DR16 sample. Each
star is color-coded by its eccentricity, with higher eccentricity stars
appearing above lower eccentricity stars. The solid line shows the separation
between halo and disc stars. The dashed shows the separation between thin and
thick disc stars. The grey contours show density in the more populated thin
and thick disc sequences.

To add a degree of realism to the _Survey_ samples we perturb all (thin disc,
thick disc, and halo) of the phase space samples using the uncertainties of
the Gaia and APOGEE DR16 data to which these samples are matched. We convert
the samples into observer coordinates: (RA, Dec., $\mu_{\mathrm{RA}}$,
$\mu_{\mathrm{Dec}}$, $d$, $v_{\mathrm{los}}$), and then add to each of the
vectors an amount sampled from a multivariate Gaussian with zero mean and
covariance matrix composed of the uncertainties and correlations drawn from
the Gaia DR2 (to stay faithful to the data used by Leung & Bovy (2019a) to
derive distances) and APOGEE data. We assert that there are no correlations
between APOGEE line-of-sight velocities and any of the Gaia data. For the
distances we use the uncertainties predicted by astroNN, and set the
correlations with other parameters to zero.

###  3.3 Sample Kinematics

For both the _Solar_ and _Survey_ set of samples we calculate kinematic
quantities as well as actions using galpy (Bovy, 2015). Actions are calculated
using the “Stäckel fudge” method of Binney (2012). In this approach our Milky
Way potential is locally approximated for each individual sample orbit as a
Stäckel potential using a focal length estimated with equation (9) of Sanders
(2012). Eccentricities are similary calculated using a variation of the
“Stäckel fudge” method as described in Mackereth & Bovy (2018). These actions
and eccentricities are calculated in the axisymmetric MWPotential2014
potential.

We consider six kinematic spaces which are commonly used throughout the
literature for both structure identification and classification, particularly
of halo populations. The first two use only velocity information, while the
latter four rely on conserved quantities calculated using the underlying
potential. For each space we provide a few recent example studies which have
made use of it, but the list is by no means exhaustive. The first is the
$v_{R}-v_{T}$ plane (e.g. Belokurov et al., 2018; Fattahi et al., 2019;
Lancaster et al., 2019; Mackereth et al., 2019a; Belokurov et al., 2020;
Feuillet et al., 2020), which is useful since disc and halo populations occupy
distinct loci defined by their mean $v_{T}$, and disc populations have much
lower tangential and radial velocity dispersions than halo populations.
Additionally the degree of radial bias in a halo population is easy to gauge
by the radial extent of its velocity ellipsoid. The second space is the Toomre
diagram (e.g. Hawkins et al., 2015; Helmi et al., 2018; Koppelman et al.,
2019b; Feuillet et al., 2020; Cordoni et al., 2020) which plots the quadrature
$\sqrt{v_{R}^{2}+v_{z}^{2}}$ of the $v_{R}$ and $v_{z}$ velocities (also known
as $v_{\mathrm{perp}}$) as a function of $v_{T}$. Here we use cylindrical
velocities, yet the Toomre diagram is also often expressed using $UVW$
velocities. In samples that span a large range in Galactic azimuth, the
$v_{R}$ and $v_{T}$ velocities project onto $U$ and $V$ in a manner that
renders the diagram less useful. The advantage of the Toomre diagram is that
it factors in $v_{z}$ in addition to $v_{R}$, further isolating halo
populations thanks to the relatively cold vertical velocity kinematics of the
thin and thick disc populations.

The next three spaces all use angular momentum about the disc rotation axis,
$L_{z}$, as the abscissa. $L_{z}$ is an effective coordinate since disc stars
tend to have characteristically large, positive values. Additionally, because
it is one of three actions, it is conserved under adiabatic conditions. The
third space is the $E-L_{z}$ plane (e.g. Helmi et al., 2018; Koppelman et al.,
2019b; Naidu et al., 2020; Feuillet et al., 2020; Cordoni et al., 2020; Horta
et al., 2021), which is perhaps the most widely used kinematic space for
structure identification in the Gaia era. This is due to the relative ease
with which these parameters are calculated, their approximate conservation
through cosmic time, the fact that the disc occupies a compact locus at high
$L_{z}$ and moderate $E$, and that the energy encodes both radius and velocity
information. We plot all energies offset by the potential of MWPotential2014
at infinity, $\Phi_{0}$, so that $E<0$ are bound orbits. The fourth space is
the root of the radial action, $\sqrt{J_{R}}$, as a function of $L_{z}$ (e.g.
Trick et al., 2019; Feuillet et al., 2020; Matsuno et al., 2021; Horta et al.,
2021). By using the radial action this space can very efficiently separate
halo stars on tangential versus radial orbits. This is because the action
labels the entire orbit while velocities, for example, will change throughout
an orbit. We use the square root of the radial action to compactify the
vertical axis for efficient visualization. The fifth space is eccentricity
($e$) as a function of $L_{z}$ (e.g. Cordoni et al., 2020). While this
specific two-dimensional space is not commonly used, eccentricity as a single
parameter is often used as a metric to separate GE from the remainder of the
stellar halo since it reflects the overall radial extent of an orbit (e.g.
Belokurov et al., 2018; Mackereth & Bovy, 2020; Naidu et al., 2020).

The final space is the “action diamond” (e.g. Vasiliev, 2019; Myeong et al.,
2019; Monty et al., 2020; Cordoni et al., 2020; Naidu et al., 2020) which in
our implementation plots $J_{z}-J_{R}$ as a function of $L_{z}$, the azimuthal
action. But both quantities are scaled by the total action defined as
$J_{\mathrm{tot}}=\lvert J_{R}\rvert+\lvert J_{z}\rvert+\lvert L_{z}\rvert$,
which causes the appearance of a characteristic diamond-shaped boundary (Note
that if the normalizing total action is
$J_{\mathrm{tot}}=\sqrt{J_{R}^{2}+J_{z}^{2}+L_{z}^{2}}$ the boundary is
circular, e.g., Naidu et al. 2020). One advantage of this space is that the
quantities are scaled by their total action, which means that disc populations
whose action budget is dominated by $L_{z}$ tend to cluster tightly in the
right corner of the bounding diamond. For the same reason halo sub-populations
which arise from, for example a single merger event, are likely to cluster at
a single location defined by the actions of the progenitor. The downside of
this space is that since all quantities are normalized by the total action,
populations which differ by an overall magnitude in their actions, such as
eccentric disc and halo populations, may appear indistinct. To combat this
effect it is often adviseable to use this space only with a sample confined to
a narrow range of radii or $L_{z}$.

Figure 3 shows these six kinematics spaces for the _Solar_ samples. The disc
samples, both thin and thick together, are shown as the black contours, while
the halo samples are shown as blue points. Each of the columns shows one of
the six different halo models. Figure 4 is the same as Figure 3, but shows the
kinematics for the _Survey_ samples instead. The disc sample contours are
constructed by binning the data with 40 bins across each dimension in each
space. The contour levels are then placed at 1, 10, and 100 samples per bin.

Figure 3: Kinematics of the six _Solar_ samples. Each of the columns shows one
of the subsets of the _Solar_ sample, which vary by the $\beta$ of the halo
component (labelled). The blue points show the halo samples while the black
contours show the disc samples, combining both thin and thick. Top row: Radial
and tangential velocities; second from top: Toomre diagram; third from top:
energy and $L_{z}$; third from bottom: $\sqrt{J_{R}}$ and $L_{z}$; second from
bottom: eccentricity and $L_{z}$; bottom row: action diamond. The units
(unlabelled) are as follows: velocities are in km s-1, actions and angular
momenta are in kpc km s-1, and energies are in km2 s-2. Figure 4: Same as
Figure 3 but for the six _Survey_ samples.

##  4 Exploring the Kinematic Properties of our Samples

In this section, we examine in detail each of the kinematic spaces we
presented in the previous section. We first provide a qualitative overview of
the kinematic spaces, especially with regards to how kinematics change with
variations in the value of $\beta$ of the stellar halo DF.

###  4.1 A Qualitative Assessment of the Kinematic Spaces

####  4.1.1 The Solar Samples

We first take a closer look at the _Solar_ sample kinematics from Figure 3.
Starting with the $v_{R}-v_{T}$ velocities, we can immediately see a strong
degree of evolution in the morphology of the halo samples as $\beta$ changes.
Specifically, as $\beta$ increases the dispersion along the $v_{R}$ axis
increases and the dispersion along the $v_{T}$ axis decreases, broadly in
accordance with the definition of $\beta$.

One subtle feature is that as $\beta$ increases the density of samples along
the $v_{T}=0$ line forms a broad plateau around $v_{R}=0$, rather than the
distinct density peak that can be seen when $\beta=0$ for example. This
behaviour is the result of a tendancy towards bimodality in the $v_{R}$
direction for high-$\beta$ samples, since stars are on more radial orbits and
therefore preferentially have non-zero radial velocities (with equal
preference to positive and negative velocities). Examining the fiducial model
highlights why it is important to keep this fact in mind. When low and
high-$\beta$ populations are shown together, it may appear to be a
superposition of two centrally concentrated ellipsoids, yet in reality the
lower-$\beta$ population provides the peak at $(v_{R},v_{T})=(0,0)$ while the
radially biased population tends to more evenly populate the low and high
radial velocity region. For much higher $\beta$ or when the density profile
allows for apocenters at much larger Galactocentric radii, one will also
expect to observe a genuine bimodality in the radial velocity direction, with
a distinct lack of stars with zero velocity.

This is noteworthy when modelling the local stellar halo velocity field using,
for example, a Gaussian mixture model. One would expect a two-component model
to fit the radially anisotropic halo better than a single component. This
assessment is supported by the results of Belokurov et al. (2018) and Fattahi
et al. (2019), who find that using two Gaussian components for the local halo
(i.e. one for the metal poor, more isotropic halo, and one for the radially
anisotropic halo) produces best-fits which have systematically negative model
residuals near $v_{r}=0$ (note they use spherical polar coordinates), and
positive model residuals towards larger radial velocities. This indicates that
the single component fit to the radially anisotropic population expects more
stars near $v_{r}=0$ and fewer at large radial velocities, in accordance with
our expectations. Other studies which have used symmetric, bimodal Gaussians
have found good fits to data (Lancaster et al., 2019; Necib et al., 2019).

The Toomre diagram shows evolution with increasing $\beta$ in a manner similar
to the $v_{R}-v_{T}$ plane. As the radial anisotropy of the halo population
increases the distribution elongates along the $v_{\mathrm{perp}}$ axis, and
narrows along the $v_{T}$ axis. The inclusion of $v_{z}$ information in
addition to $v_{R}$ in the Toomre diagram does tend to separate the halo
samples from the disc samples somewhat more than in the $v_{R}-v_{T}$ plane.
But the effect is not as pronounced as might be expected, likely because many
disc samples have roughly comparable vertical and radial velocities, and so
their extent along the perpendicular velocity axis also increases. In both of
these first two kinematic spaces, the disc occupies similar loci at $v_{T}\sim
200$ km s-1 with zero net radial or vertical velocity. It is important to
note, however, that even in this population of disc samples located at the
solar position there are a marginal yet noticable number belonging to the
thick disc with near-0 $v_{T}$. In the Toomre diagram even though the halo
samples are more distinct, the maximum perpendicular velocities of the hottest
thick disc samples approaches 200 km s-1.

The $E-L_{z}$ plane in the third row of Figure 3 exhibits the now famous
central enhancement of samples with $L_{z}\approx 0$ for large $\beta$. It is
this feature which is most characteristic of the GE population as first shown
by Helmi et al. (2018). One of the reasons this space has been studied so
frequently in the Gaia era is that both $E$ and $L_{z}$ are approximately
conserved integrals of motion, and therefore the appearance of stellar
(sub)structures in this space often point to a common ancestry for the
constituent stars. Another key feature of this space which we will explore in
more depth later is the distinct parabolic shape of the overall distribution.
This boundary is set by the fact that for a given energy and radius the
magnitude of the velocity is fixed. This results in a bounding value of
$L_{z}$ which occurs when the velocity is wholly oriented in the direction of
$v_{T}$. The equation of the bounding curve is:
$E=L_{z}^{2}/(2R^{2})+\Phi(R)$. Because the samples in the _Solar_ set are all
at a fixed radius the boundary is obvious in Figure 3. There is also a clear
gradient in energy, with the density of samples decreasing as the energy
increases. This is intrinsic to the energy-dependence of spherical DFs (see
chapter 4.3 in Binney & Tremaine, 2008) and is therefore observed even at a
fixed radius. As the radius changes, the energy of the apex of the bounding
parabola also changes. This emphasizes that the relationship between radius
and energy is complicated, a feature which we will discuss more later.

The radial action, $J_{R}$, is perhaps the most natural variable for isolating
stars from populations with high radial anisotropy, because it labels the
entire underlying orbit rather than being sensitive to the phase of the orbit
such as velocities are. We see in Figure 3 that there is an approximately
inverse relationship between the magnitude of $L_{z}$ and the radial action.
This occurs because as angular momentum increases at a fixed position, less of
the velocity budget is necessarily contained in radial motions. The
relationship is not perfect, however, for two reasons. First, this kinematic
space does not account for vertical motions, and second there is a range of
energies in the _Solar_ samples as we have just discussed. As $\beta$
increases a distinct cusp forms around $L_{z}=0$ spanning a range of
$\sqrt{J_{R}}$. The fact that this cusp has a distinct minima (around
$\sqrt{J_{R}}\approx 25$
$\mathrm{kpc}^{1/2}~{}\mathrm{km}^{1/2}~{}\mathrm{s}^{-1/2}$ in Figure 3)
reflects the fact that at the solar circle, the only orbits that exist with
$L_{z}=0$ must either have a certain degree of radial motion, or otherwise
must be on a circular orbit exactly over the Galactic pole. The cusp spans a
range of $\sqrt{J_{R}}$ for the same reason that the inverse relationship
between $\sqrt{J_{R}}$ and $L_{z}$ is imperfect at non-zero $L_{z}$: there is
a range of energies in the samples and vertical motions are not accounted for.

The relationship between eccentricity and $L_{z}$ is similar in spirit to the
relationship between $\sqrt{J_{R}}$ and $L_{z}$. As $L_{z}$ increases at a
fixed location, eccentricity must necessarily decrease because orbits become
less radial. This lends a characteristic v-shape to the samples in this space.
Unlike $\sqrt{J_{R}}$, eccentricity is a “scaled” parameter, it depends only
on the relative values of the orbit’s pericenter and apocenter. Therefore, we
see that disc samples appear less distinct from halo samples than in other
kinematic spaces. From an intuitive perspective, a halo star with pericenter
at the solar circle could have the same eccentricity as a disc star with
apocenter at the solar circle, yet these two stars would have very distinct
radial actions. A second curiosity of eccentricity is the fact that as orbits
become more radial, their eccentricity approaches 1 asymptotically slowly.
This means that as $\beta$ increases samples bunch up close to $L_{z}=0$ and
$e=1$. This is markedly different behaviour than most other kinematic spaces
we have presented so far, where the kinematic properties of radially-biased
samples have a distinct magnitude. Here, the only distinct property of the
radially biased population is its sheer density at high eccentricity. We will
discuss how this impacts the use of eccentricity as a parameter to select
high-$\beta$ populations later.

The action diamond is perhaps the least used among the kinematic spaces shown
here, with a similar implementation first being used by Vasiliev (2019) (to
our knowledge, although they recognize that the space is analogous to
visualisations of three-dimensional action spaces presented by Binney &
Tremaine, 2008). Its most common use thus far has been to study globular
cluster kinematics and their potential association with larger halo structures
(Vasiliev, 2019; Myeong et al., 2019). But it still offers useful information
for the types of aggregate stellar populations we study here. Understanding
the action diamond is intuitive. The right and left points of the diamond
($\lvert L_{z}\rvert=J_{\mathrm{tot}}$) are in-plane prograde and in-plane
retrograde orbits respectively. The top point ($J_{z}=J_{\mathrm{tot}}$) is a
polar orbit, and the bottom point ($J_{R}=J_{\mathrm{tot}}$) is a radial
orbit. The bottom-right and bottom-left edges ($J_{\mathrm{tot}}$ wholly
contained in $L_{z}$ and $J_{R}$) are prograde and retrograde in-plane orbits.
The top-right and top-left edges ($J_{\mathrm{tot}}$ wholly contained in
$L_{z}$ and $J_{z}$) are prograde and retrograde out-of-plane orbits.

One of the most noticeable features of the action diamond is that both disc
and halo populations tend to be biased towards the bottom half of the diamond,
which reflects the fact that both populations tend to have modestly hotter
radial versus vertical kinematics. Indeed, we see what is perhaps the most
subtle transformation in morphology from the low to high-$\beta$ halo
population, where the samples appear to concentrate more towards the bottom
point of the diamond rather than the left and right points. This results from
the dominance of the radial action over the vertical action in the
high-$\beta$ sample. The disc also hugs the bottom-right edge of the diamond,
but the action budget is dominated by $L_{z}$, meaning it is mainly confined
to the right point of the diamond. Similarly to eccentricity this is also a
“scaled” space, so two orbits with different absolute kinematic properties may
occupy the same point. So, while the strength of the action diamond lies in
the fact that distinct types of orbits occupy different loci in the space, the
weakness is that confusion between two populations that may have different
absolute kinematics, but similar types of orbits is likely. This will become
more evident when examining the _Survey_ sample set.

####  4.1.2 The Survey Samples

The _Survey_ samples differ from the _Solar_ samples in three key respects
with regards to the six kinematic spaces that we consider. First, the sample
positions span a range of Galactocentric radii, and so in some sense it is
useful to consider the kinematic profiles as a superposition of many of the
individual samples seen in Figure 3, but determined at different radii.
Second, the samples are perturbed by Gaia and APOGEE errors. Third, and
perhaps most important, is that the positions underlying the samples are
subject to the complicated APOGEE DR16 selection function. It is through these
three lenses that we will assess the kinematics of the _Survey_ sample set in
Figure 4.

First, one obvious change to all of the kinematic spaces is the increased
range over which the disc populations are observed, while in the _Solar_
samples the disc populations tended to be confined to a narrow region. This is
due mostly to the first two factors mentioned above: the increased radial
range over which the DFs are sampled and the addition of realistic
uncertainties. Furthermore the discrepancy in the number of disc samples
versus halo samples (roughly 66:1) exacerbates the degree to which the disc
populations extend into the stellar halo kinematic space with meaningful
density.

In the $v_{R}-v_{T}$ plane, the overall morphology of the halo samples stays
the same, with bimodality along the $v_{R}$ axis increasing as $\beta$
increases. The impact of the realism steps on the data is clear, particularly
in how the _Survey_ samples do not cluster so tightly on the $v_{T}=0$ line,
especially in the samples with higher $\beta$, as the _Solar_ samples do.
Overall our fiducial model looks realistic in its representation of the nearby
Milky Way $v_{R}-v_{T}$ plane (c.f. figures 2 and 1 in Belokurov et al., 2018;
Fattahi et al., 2019, respectively, although note they use spherical
coordinates), a testament both to the validity of our models and the steps we
take to add realism to our DF samples. The Toomre diagram changes in much the
same way as the $v_{R}-v_{T}$ plane: the overall morphology is the same, but
the addition of uncertainties has softened the cusp in the halo samples along
the $v_{T}=0$ line. In both these planes, the disc population now extends down
to near $v_{T}=0$, which reflects the addition of samples from the inner
galaxy where the angular momentum is lower and the radial velocity dispersion
is higher. The extent of the disc in the $v_{R}$ and perpendicular velocity
dimensions does not increase as drastically, because only the increased radial
and vertical velocity dispersions in the inner Galaxy serve to modify its
extent.

From a qualitative standpoint, the $E-L_{z}$ plane changes in much the same
way as the velocity planes. The cuspy overdensity of the high-$\beta$
populations is softened, and the hottest thick disc stars begin to impinge on
this characteristic region, yet leave a discernable overdensity of halo
samples visible just beyond their boundary. The $E-L_{z}$ plane also shows a
distinct change in morphology of the bounding region of the kinematic space,
which we mentioned in the previous section is due to there being a fixed
velocity budget at a given energy and radius. For the _Survey_ sample, where
stars at many radii are included, the bounding region is a superposition of
many of the aforementioned parabolic curves. These bounding regions tend to
get wider as radius increases since the potential scales inversely with
distance from the Galactic center. For inner Galaxy samples the bounding
parabola is narrow with minima at very negative energies. Conversely, for
outer Galaxy stars, the bounding parabola is wide with a minima at higher
energy. The superposition of these bounding curves results therefore in a
funnel shape which narrows as the energy is decreased. An example of this
morphology is shown in Figure 5 where the in-plane bounding parabolas for
MWPotential2014 corresponding to radii 2-15 kpc are shown.

Figure 5: The $E-L_{z}$ plane shown with bounding parabolas corresponding to
Galactocentric radii 2 (lowest curve) to 15 (highest curve) kpc. The curves
are calculated in-plane ($z=0$), in MWPotential2014.

A particularly interesting feature of the _Survey_ sample $E-L_{z}$ plane is
that there appears to be a region of maximal density which begins at
$E\sim-1.3\times 10^{5}$ km2 s-2 and follows an approximate parabolic shape
towards larger energies as $\lvert L_{z}\rvert$ increases. These are most
obvious in the $\beta=0$ and $0.5$ models, but are also noticable at higher
$\beta$. Below this energy there are halo samples, yet they are observed in
much lower numbers than at this overdensity. The interpretation of this
overdensity is that the majority of the samples in the _Survey_ set lie near
the Sun, and therefore they tend to be roughly confined by the bounding
parabola defined at that location. The occurence of over- and underdensities
in $E-L_{z}$ space depending on the underlying radial distribution of samples
is something we will explore in the next subsection.

The $\sqrt{J_{R}}-L_{z}$ plane no longer has any indications of the inverse
relationship between $J_{R}$ and $L_{z}$ observed in the _Solar_ samples and
discussed above. Likely the combination of the radial range as well as
realistic uncertainties causes all traces of the relationship to disappear.
Instead now the halo samples manifest as a broad cloud that is slightly biased
either towards $L_{z}$ or $\sqrt{J_{R}}$ depending on the value of $\beta$.
The halo samples of the most radially biased models still do exhibit a diffuse
cusp along the $L_{z}$ = 0 line. The differences in the $e-L_{z}$ plane tell a
similar story. There is no longer as sharp an inverse relationship between
eccentricity and $L_{z}$ as was observed in the _Solar_ samples, yet the
relationship is at least somewhat apparent in contrast with the
$\sqrt{J_{R}}-L_{z}$ plane. The higher $\beta$ models still tend to simply
cluster more tightly around $L_{z}=0$ and $e=1$, similar to the _Solar_
samples. Notably, the disc samples extend all the way to the high-$\beta$
overdensity, and in general completely dominate the prograde half of the
plane.

In the action diamond, the disc populations transition from occupying a narrow
strip near the prograde corner to dominating the lower right quadrant of the
diamond. This large change in area of the disc populations demonstrates a
weakness in the action diamond (indeed one shared by eccentricity), which is
that “scaled” kinematic quantities often confuse orbits that only broadly
share certain properties. In this case, when the vertical and radial actions
of the hottest thick disc stars become of order comparable with their $L_{z}$
they migrate towards the center and bottom corner of the diamond overlapping
with the halo samples, even though the magnitude of their actions can be quite
different from those of stars in the halo. The halo samples have the same
qualitative trends as are seen in the _Solar_ samples. As $\beta$ increases,
the samples concentrate near the bottom of the figure. Additionally of note is
that very few ergodic samples occupy this region of the diagram, they are more
uniformly spread throughout the space and even biased towards the high-$L_{z}$
corners. This means the bottom corner is likely a good region to isolate
radially biased halo populations from ergodic populations.

###  4.2 The Impact of the Underlying Radial Distribution of Samples on the
$E-L_{z}$ plane

In the previous section, we noted that the _Survey_ samples exhibit an
interesting change in morphology in the $E-L_{z}$ plane when compared with the
_Solar_ samples. In the $E-L_{z}$ plane of the latter (Figure 3) the samples
were confined by a bounding parabola defined by the escape velocity as a
function of changing energy at a fixed Galactocentric radius (this being the
location of the Sun in the _Solar_ sample). Sample density decreases as energy
increases moving away from the bounding parabola. The radially biased samples
have a higher energy cusp around $L_{z}=0$, leading to different overall
morphology but the same trend still holds. This phenomenon lead us to
interpret the apparent overdensity in the $E-L_{z}$ plane of the low-$\beta$
_Survey_ samples around $E\sim-1.3\times 10^{5}$ km2 s-2 as being caused by
the concentration of the samples near the Sun.

While this observation is not particularly unexpected, it suggests a question:
what would happen if the intrinsic radial distribution of samples (i.e. the
APOGEE DR16 data which sets the approximate positions of our samples) took a
more complicated form? We can investigate this using APOGEE DR16 by augmenting
the _Survey_ sample with additional samples directed towards the Galactic
bulge. Recall that we do not consider any stars for inclusion in the _Survey_
sample if they have $-20\degr<\ell<20\degr$ and $|b|<20\degr$. Here, just to
investigate the radial distribution of samples, we soften that requirement and
allow stars lying in this angular range with Galactocentric radius greater
than 3 kpc to be included into the _Survey_ sample set. Note that this
distance from the Galactic center is approximately equal to the distance from
the Galactic center enforced by the $20\degr$ angular cut we use to define the
_Survey_ sample - i.e. $\tan(20\degr)\times 8.178~{}\mathrm{kpc}\approx
3~{}\mathrm{kpc}$. This radius is also the radius at which the density of the
bulge component of MWPotential2014 drops below the density of our stellar halo
model when it is normalized to have the mass found by Mackereth & Bovy (2020).
This is all to say that we are confident we do not introduce many bulge stars
into the sample by including these new data. Separating halo stars from disc
stars, sampling kinematics from the DF, and the perturbation by uncertainties
are all done in the same manner as described in § 3. This results in 725 new
halo samples.

The top panel of Figure 6 shows the radial distribution of the _Survey_
samples, the additional samples in the direction of the Galactic center, and
their combination. Note that a number of the additional Galactic center
samples appear at radii less than 3 kpc. This can be attributed to the
perturbation by the APOGEE and Gaia uncertainties. The overall distribution of
radii is clearly bimodal, with peaks around 3.5 and 8 kpc, and a decrease at 6
kpc. The inner peak is driven mostly by the additional samples from the inner
Galaxy pointings, but hints of the overdensity are clearly present in the
_Survey_ sample. The occurence of this bimodality is undoubtedly caused by the
APOGEE selection function, and specifically the location of the survey
pointings.

Figure 6: Top: Histogram of Galactocentric radii for the _Survey_ sample
(blue), the additional samples towards the Galactic bulge (red), and their
combination (black). Bottom left: $E-L_{z}$ plane for _Survey_ sample
augmented by the additional sources towards the Galactic bulge. Bottom right:
Histogram of energies for the same sample. Three grey dashed lines in the top
panel mark the locations of over and underdensities in the radial distribution
at 3.5, 6, and 8 kpc. Three black lines in the bottom panel show sections of
the bounding parabolas for those corresponding radii, roughly marking the
boundaries between over and underdensities in the $E-L_{z}$ plane.

The bottom-left panel of Figure 6 shows the $E-L_{z}$ plane for the fiducial
model of the _Survey_ sample plus the additional samples towards the Galactic
center, 3523 samples in total (halo samples only, no disc). It is immediately
clear that there is an underdensity between $-1.25$ and $-1.5\times 10^{5}$
km2 s-2. This is confirmed by examining the bottom right panel, a histogram of
energies, which highlights the overdensity at the aforementioned location. It
does not perfectly highlight the overdensity because, as we have established,
the density fluctuations occur on roughly parabolic contours, but nontheless
its signature is there. Overlaid on the figure are three black lines that are
sections of the bounding parabolas corresponding to the three radii marked
with grey dashed lines in the top panel: 3.5, 6, and 8 kpc. The bottom and top
curves (3.5 and 8 kpc) cradle the two overdensities, while the middle curve at
6 kpc cradles the underdensity. This clearly illustrates that the uneven
radial distribution of samples is clearly linked to the appearance of
substructure in the $E-L_{z}$ plane.

##  5 The Purity of Populations of Interest

One of the most useful questions we might ask is whether we can use our
samples to identify which kinematic spaces are best for identifying certain
populations that we are interested in separating from a background of
contaminants. For the purposes of this work we focus on separating the
high-$\beta$ populations from low-$\beta$ and disc populations.

###  5.1 Separating High-$\beta$ Populations from Low-$\beta$ and Disc
Populations

The separation of high-$\beta$ populations from low-$\beta$ populations is
currently of special importance, because the two dominant inner-Milky Way halo
populations are the in-situ low-$\beta$ component and the high-$\beta$ GE
component (Belokurov et al., 2018; Helmi et al., 2018; Haywood et al., 2018).
While these two populations do have somewhat unique chemistry (Haywood et al.,
2018), it is not so distinct that that abundances are an efficient means of
separating them, such as is the case for separating the disc populations both
from one another and from the halo. Building on this sentiment, we do not even
know a priori the detailed chemical profile of the in-situ and GE components
so that we might try and separate them using their abundances. Their
individual chemistries can only be effectively learned if the populations can
be separated by some other means. Additionally, detailed mass modelling in the
halo requires some means of tagging different populations of interest. For
example Mackereth & Bovy (2020) divide the halo into mono-abundance
populations. They use a simple eccentricity cut to separate ergodic and
high-$\beta$ populations. Here we will delve deeper into the methods for
kinematically separating these two major halo populations.

We use our fiducial halo model to investigate the purity of ergodic and
high-$\beta$ populations; recall that the fiducial halo model consists of a
50:50 mixture of $\beta=0.5$ and 0.9 populations. Figures 7 and 8 show the
purity of the high-$\beta$ component of the fiducial model with respect to the
low-$\beta$ component for the _Solar_ and _Survey_ sample sets respectively.
Purity is calculated as
$N_{\mathrm{\beta=0.9}}/(N_{\mathrm{\beta=0.9}}+N_{\mathrm{\beta=0.5}})$, and
the data is binned with 40 bins across each dimension of each kinematic space.
Overlaid on each kinematic space is the lowest contour of the disc samples
from Figures 3 and 4 to indicate where the disc may contribute contaminants to
the purity measurement.

Figure 7: Purity of high-$\beta$ halo samples with respect to low-$\beta$
samples for the _Solar_ sample set. The kinematic spaces and units are as
described in Figure 3. In each of the panels the dashed black ellipse(s) show
the selection ellipses for the high-$\beta$ sample. The solid black contour is
the lowest contour of disc stars from Figure 3. Figure 8: Same as Figure 7 but
for the _Survey_ sample set. The solid black contour is the lowest contour of
disc stars from Figure 4.

Qualitatively, the purity of high-$\beta$ samples follows intuitively from
Figures 3 and 4. In the $v_{R}-v_{T}$ and Toomre diagrams the purity is
highest along the $v_{T}=0$ axis and at large values of $v_{R}$ and the
perpendicular velocity respectively. In the $E-L_{z}$ plane, we see that the
purity is highest when $L_{z}$ is zero and at modest to high energies. The
$\sqrt{J_{R}}-L_{z}$ plane behaves similarly, with a cloud of high-purity at
large values of the radial action near $L_{z}=0$. Eccentricity is peculiar, as
we have already noted in that being a “scaled” parameter the radially biased
populations do not occupy a distinct region of the space. Rather they tend to
cluster with greater density near $L_{z}=0$ and $e=1$. This behaviour is
necessary to keep in mind when evaluating this space on the merits of
high-$\beta$, because one might assume that the best place to isolate radially
biased populations is near $e=1$ but at higher $L_{z}$. In reality, the bulk
of the sample lies at $L_{z}=0$ but simply exists at lower purity because the
ergodic sample also populates that region. The action diamond exhibits the
highest purity at the bottom corner where the radial action dominates the
action budget. In particular in the _Solar_ sample, the region of high-purity
is only seen in the few pixels around the apex of the diamond.

Table 2: High-$\beta$ selection ellipse parameters for the _Solar_ and _Survey_ samples in six kinematic spaces. Ellipses correspond to those in Figures 7 and 8. See the caption of Figure 3 for unit conventions. From left to right the columns are: kinematic space, ellipse center point, ellipse semi-major axes, high-$\beta$ completeness, high-$\beta$ purity for halo samples only, and high-$\beta$ purity including disc samples. Space | Center [x,y] | SMA [x,y] | C | Phalo | Pdisc
---|---|---|---|---|---
_Solar_ sample |
$v_{R}-v_{T}$ | $[\pm 260,0]$ | [140,20] | 0.28 | 0.80 | -
Toomre | $[0,200]$ | [20,200] | 0.72 | 0.74 | -
$E-L_{z}$ | $[0,-0.9]$ | [200,0.45] | 0.71 | 0.75 | -
$\sqrt{J_{R}}-L_{z}$ | $[0,45]$ | [250,20] | 0.69 | 0.78 | -
$e-L_{z}$ | $[0,1]$ | [500,0.05] | 0.70 | 0.77 | -
Action Diamond | $[0,-1]$ | [0.1,0.1] | 0.62 | 0.86 | -
_Survey_ sample |
$v_{R}-v_{T}$ | $[\pm 280,0]$ | [100,40] | 0.10 | 0.73 | 0.71
Toomre | $[0,280]$ | [35,100] | 0.12 | 0.69 | 0.68
$E-L_{z}$ | $[0,-1]$ | [300,0.5] | 0.65 | 0.65 | 0.60
$\sqrt{J_{R}}-L_{z}$ | $[0,45]$ | [300,15] | 0.18 | 0.76 | 0.75
$e-L_{z}$ | $[0,1]$ | [500,0.05] | 0.48 | 0.74 | 0.72
Action Diamond | $[0,-1]$ | [0.08,0.3] | 0.23 | 0.82 | 0.80

When evaluating how each kinematic space performs in separating populations,
it is important to keep in mind the overall density of halo stars. Many
regions of each space exhibit spuriously high purity where there are only a
few halo samples. This effect is more pronounced in the _Survey_ samples than
in the _Solar_ samples, because of the decreased numbers (2,798 versus
$10^{4}$). With this in mind, we have selected elliptical regions in each
space where both the purity is large and the density of high-$\beta$ halo
samples is large. These are shown as the dashed black lines in Figures 7 and
8. The parameters of the selection ellipses are given in Table 2. We choose
these simple shapes to represent the high-$\beta$ regions because one goal in
this analysis is to emphasize applicability to real data. Shapes which are too
specific to the populations here may provide errant results when applied to
other data where the kinematic properties may be subtly different (e.g., if
either the spatial selection or the potential used to derive kinematics
changes). Instead, we pick regions which would be mostly insulated from the
peculiarities of a specific data set. In each of the six kinematic spaces at
least one of the dimensions probes the circularity of the orbit in the disc
plane. Following this, each of the selection regions primarily picks out
samples with $L_{z}$ or $v_{T}$ approximately equal to zero (the bottom corner
of the action diamond necessarily has $L_{z}=0$). The selections then pick out
extreme values of the other quantity, which generally trends with the degree
of radial bias in a population (specifically $v_{R}$, perpendicular velocity,
$E$, and $\sqrt{J_{R}}$). The selection regions of the “scaled” spaces tend to
be concentrated at the extrema of the scaled quantity (eccentricity and
$(J_{z}-J_{R})/J_{\mathrm{tot}}$). The reason for mentioning the overall
trends in our choice of selection regions is that these broad properties
should suffice to pick out radially biased samples in any data set, not just
specifically our DF models. This is important if, for example, we are to be
able to confidently apply these selections to real APOGEE DR16 data, something
we plan to do in the near future.

It is still interesting to ask, however, what is the tradeoff between
increasing or decreasing the size of the selection ellipses? Intuitively, a
larger selection boundary would include more high-$\beta$ samples, but would
also include more low-$\beta$ samples which would decrease purity. To explore
this, we make purity-completeness relations for each of the kinematic spaces
for both the _Solar_ and _Survey_ sample sets. To do this we consider two sets
of scaling factors: $(x_{1},x_{2})$ where each factor takes the form
$x_{i}=2^{y}$. For each of the two scaling factors we consider 20 values of
$y$ uniformly spaced on the interval $[-2,2]$, giving a grid of 400 total
combinations. For each of these 400 unique pairs of factors we multiply them
to the semi-major axes in Table 2 to generate a new selection ellipse. We then
calculate the high-$\beta$ purity and completeness according to this new
selection ellipse. To generate an average purity-completeness relationship we
bin the 400 purity measurements in completeness and calculate the mean purity
in each bin. The results are shown for the _Solar_ and _Survey_ sample sets in
Figures 9 and 10 respectively. Each of the panels highlights one kinematic
space (labelled), the curves for the other kinematic spaces are shown faintly
behind to emphasize differences. In the panels of Figure 10 there is a solid
and dashed line that show the purity calculated only considering halo samples
(solid) and the purity calculated including both halo and disc samples
(dashed). Additionally, in each panel there is a colored circled which shows
the purity and completeness of the default boundaries shown in Figures 7 and
8.

Figure 9: Purity-completeness relation for the fiducial _Solar_ sample for
each kinematic space (labelled). The bold line in each panel shows the purity
calculated only considering halo samples, the solid lines from all other
panels are shown faintly in the for context. The filled circle in each panel
shows the purity and completeness of the default selections from Figure 7.
Figure 10: Same as Figure 9 but for the fiducial _Survey_ sample.
Additionally, the bold dashed line in each panel shows the purity calculated
including disc samples. Also see Figure 8.

For the _Solar_ samples, the purity-completeness relations are similar for
quantities which are similar to one another. $v_{R}-v_{T}$ and the Toomre
diagram are similar, $E-L_{z}$ and $\sqrt{J_{R}}-L_{z}$ are similar, and
$e-L_{z}$ and the Action Diamond are similar. The Action Diamond appears to be
the superior space, achieving a higher purity at a given completeness than any
other kinematic space. At large completeness the other kinematic spaces are
more-or-less tied with one another. The apparant bimodalities (i.e. the peaks
near completeness of 0.6) in the purity-completeness relations for the
absolute kinematic spaces is just an artifact of how we have visualized the
results of this analysis. In general, the linear trend at higher completeness
in these spaces captures the variation in the rotation-sensitive quantity
($v_{T}$ and $L_{z}$) when the anisotropy-sensitive quantity ($v_{R}$,
$v_{\mathrm{perp}}$, $E$, and $\sqrt{J_{R}}$ is scaled to a large range. At
lower completeness the same trend holds when the radial quantity is scaled to
small values. These two trends emerge because of the logarithmic manner in
which we generate the scaling factors. A general trend for these four spaces
is as follows: at a given ellipse size for the anisotropy-sensitive quantity,
a narrower selection for the rotation-sensitive quantity produces higher
purity and lower completeness. For the Action Diamond and $e-L_{z}$ the
relations are simple because the high-$\beta$ samples occupy a single locus in
these spaces. In general the only takeaway from this analysis applied to the
_Solar_ samples is that the action diamond appears to outperform the other
spaces. We are wary of making any other statements because the simplicity of
our _Solar_ samples clearly incites slightly pathological results in this
specific analysis.

For the _Survey_ samples the differences between the kinematic spaces are more
interesting. The action diamond and $e-L_{z}$ plane have very similar purity-
completeness relations, with the action diamond having lower completeness and
higher purity, and the $e-L_{z}$ plane having higher completeness and lower
purity. Both of these spaces have little contamination from disc samples
(solid and dashed lines are similar). These two scaled spaces appear superior
to the other four spaces under consideration. A close third though appears to
be the $\sqrt{J_{R}}-L_{z}$ plane, which has slightly lower purity, with the
other three spaces at the lowest purities. The shapes of the purity-
completeness relations for the velocity planes appears to be slightly peaked
at modest completeness. These modest peaks represent selection ellipses that
span the whole range of $v_{R}$ and $v_{\mathrm{perp}}$ respectively. With
these selections you would greatly increase the completeness while keeping the
purity relatively constant, however the purity including disc stars drops
drastically. Overall the $v_{R}-v_{T}$ plane achieves slightly higher purity
than the Toomre diagram, likely confirming earlier suspicions that the
inclusion of $v_{z}$ in the Toomre diagram’s $v_{\mathrm{perp}}$ actually
serves to confuse between low- and high-$\beta$ samples. Finally, the
$E-L_{z}$ plane achieves purities comparable to the velocity planes, but it
does so with quite high completeness.

For both the _Solar_ and _Survey_ samples we also consider composite
selections combining the selection ellipses of multiple kinematic spaces.
Intuitively, one might expect that the best composite selection would be to
combine a scaled space, such as the action diamond or $e-L_{z}$ with
$\sqrt{J_{R}}-L_{z}$ or $v_{R}-v_{T}$. The scaled space selects for stars on
very radial orbits, while the other space selects for stars with radial
motions of a large magnitude. We find only very minor improvements by
combining the action diamond or $e-L_{z}$ selections with the
$\sqrt{J_{R}}-L_{z}$ or $v_{R}-v_{T}$ selections for the _Survey_ samples. The
improvement in purity is at most 2 per cent (and sometimes purity decreases),
while it tends to cut completeness by a factor of two or more. For example,
when combining the action diamond selection with the $\sqrt{J_{R}}-L_{z}$
selection purity actually decreases by 0.01 to 0.81 (compared to just the
action diamond), and completeness falls to 0.1 from 0.23. One positive aspect
of combining selections is it does tend to erase any disparity between the
purity with and without disc samples. Nonetheless this is a major hit to
completeness for a very minor benefit in terms of purity. The story is the
same for the _Solar_ samples. With all of this information in mind, the best
composite selection (for the _Survey_ samples) appears to be the union of the
action diamond and $e-L_{z}$ selections. When compared to just the action
diamond selection, the purity without (with) disc contaminants is 0.82 (0.81)
and the completeness only falls by 0.01 to 0.22. Although this selection is
only marginally better than just using the action diamond alone. If maximum
purity is of interest, then the action diamond plus the $v_{R}-v_{T}$
selection appears best, where the purity without (with) disc contaminants is
0.84 (0.83) however the completeness is only 0.06.

This analysis does have a few restrictions. First, we impose the use of simple
ellipses because, as we argued above, they are simple selection shapes that
can be confidently applied to real data. Second, we fix the location of the
ellipses in each kinematic space. We also explored an alternative method for
gauging purity-completeness relations for the _Solar_ and _Survey_ sample sets
which does not assume a shape for the selection region, and instead works on a
pixel-by-pixel basis. The advantage of this approach is that it examines
purity and completeness in an unbiased manner. The weakness is that any
selection boundaries based on this individual pixel approach are highly
sensitive to the exact kinematic configuration of these data. With this in
mind we would not be confident using this method to define selection
boundaries for extracting high-$\beta$ stars from real data. The ellipse based
method, while biased, is sufficiently general that it should work well on real
data.

###  5.2 Eccentricity as a Single Metric for Separating the High-$\beta$ and
Ergodic Halo

In the analysis presented in the previous section, we found that good purity
of high-$\beta$ halo samples with respect to low-$\beta$ halo samples and disc
samples can be achieved using eccentricity only if a very narrow region around
$e=1$ and $L_{z}=0$ is selected. Outside of this region confusion between
high-$\beta$ samples and these contaminants is high. The reasons for needing
this very specific selection criterion is as follows: first, that eccentricity
is a “scaled” quantity, meaning two physically distinct orbits that share a
common ratio of apocenter to pericenter have the same eccentricity; second,
that as orbits become increasingly radial, eccentricity approaches unity
asymptotically. Indeed, examining the $e-L_{z}$ planes in Figures 3 and 4
shows that a reasonable fraction of the sample in the ergodic halo models lies
at eccentricity near 1.

In the literature, eccentricity is commonly used as a metric to both describe
(Myeong et al., 2019; Mackereth et al., 2019a), and separate (Mackereth &
Bovy, 2020; Naidu et al., 2020) low- and high-$\beta$ halo populations.
Oftentimes eccentricity is used by itself, as a single parameter without the
help of $L_{z}$ or another parameter that selects for tangential orbits. In
order to investigate how well these sorts of selection criteria work, we plot
the eccentricity distribution function for the halo samples of each of the six
models in the _Survey_ sample set in the top panel of Figure 11. We see, as
expected, that the peak of the distribution moves towards larger
eccentricities as $\beta$ increases from $-0.5$ to $0.9$. Interestingly, we
see that the ergodic and $\beta=-0.5$ models have a non-negligible number of
samples with very high eccentricity. The ergodic model in particular has a
nearly uniform density of samples between eccentricities 0.4 and 1, with a
minor apparent peak around eccentricity of 0.8. This suggests that while high
eccentricity cuts do a good job of preferentially picking out high-$\beta$
stars, there is a significant amount of contamination by the well-populated
high-eccentricity tail of the low-$\beta$ populations.

To quantify this confusion, we plot the completeness and purity of the
high-$\beta$ subsets with respect to the low-$\beta$ subset of the fiducial
_Survey_ sample in the bottom panel of Figure 11. These quantities are
expressed as a function of a cutoff eccentricity, below which we do not
consider any samples. The curves are constructed in this way: For each cuttoff
eccentricity we take the samples of the fiducial model which have greater
eccentricity and calculate purity as
$N_{\mathrm{\beta=0.9}}/(N_{\mathrm{\beta=0.9}}+N_{\mathrm{\beta=0.5}})$ and
completeness as $N_{\mathrm{\beta=0.9}}/N_{\mathrm{\beta=0.9,tot}}$ where
$N_{\mathrm{\beta=0.9,tot}}=$ is 50 per cent of the total number of halo
samples in the fiducial model, as per its definition. For smaller eccentricity
cuttoffs $\lesssim 0.8$ the purity and completeness vary inversely, but as the
cuttoff increases the completeness drops sharply while the purity only
continues to rise in a roughly linear fashion.

Figure 11: Top: eccentricity distribution for the halo samples of each of the
six models in the _Survey_ sample set. The dashed line marks eccentricity=0.7,
a traditional boundary between the ergodic and high-$\beta$ halo populations.
Bottom: purity and completeness of the high-$\beta$ samples with respect to
the ergodic samples in the fiducial halo model of the _Survey_ sample set. The
quantities are expressed as a function of the cutoff boundary in eccentricity.

Notably, when the eccentricity cuttoff is 0.7 the purity of the high-$\beta$
sample is only 0.59, with the completeness being 0.9. This is a remarkable
amount of confusion between these two samples. Recall that the high-$\beta$
purity of the entire fiducial sample is 0.5 by definition, and an eccentricity
threshold of 0.7 does little to improve that. When the eccentricity cuttoff is
increased to 0.9, the purity rises to 0.68 and the completeness is reduced to
0.66. Even at the largest value we consider for the eccentricity cuttoff,
0.99, the purity only reaches 0.8 (roughly in accordance with our variation of
ellipses analysis in the previous section), and the completeness plummets to
only 0.17. These results clearly indicate that care need be taken when
separating high-$\beta$ and ergodic populations on the basis of eccentricity
alone.

##  6 Discussion

The data and results we present in this paper are clearly useful to discuss in
the context of a large fraction of the literature on Milky Way stellar
dynamics done in the Gaia era, specifically those which focus on the stellar
halo. Here we focus on discussing the takeaways we deem are most salient, and
how they relate to the wide body of literature.

###  6.1 The Validity and Applicability of our Models

First, we must assess how good the models we have presented are at
representing the Milky Way halo as we understand it to be today. From a
theoretical perspective, our models are clearly a simplification. For example,
it is now well-known that the Milky Way stellar halo is effectively described
by a metallicity-dependent triaxial power law density profile (e.g., Mackereth
& Bovy, 2020). Our models at least have a power law index and cuttoff radius
which match observations. But more complicated models are limited by the
availability of usable DFs. Additionally, the Milky Way stellar halo is
unlikely to currently be in a completely relaxed state, with perturbations
from nearby massive dwarf galaxies in various stages of merging (Garrow et
al., 2020), and specifically the Large Magellanic Cloud (Erkal et al., 2019).
Fortunately for us these effects are greatest in the outer halo, with
dynamical times in the inner halo where our sample is concentrated being short
enough such that these perturbers do not meaningfully impact the equilibria.

A better method of addressing the validity of our models is to take an
empirical approach and compare the kinematics of our samples with real
observed stellar kinematics. For example, when considering the fiducial
_Survey_ model, the $v_{R}-v_{T}$ planes of Belokurov et al. (2018); Fattahi
et al. (2019); Naidu et al. (2020), the Toomre diagram of Helmi et al. (2018);
Koppelman et al. (2019b), the $E-L_{z}$ plane of Helmi et al. (2018);
Koppelman et al. (2019b); Naidu et al. (2020), and the action diamond of
(normalized to be a circle, Naidu et al., 2020) all show similarities with the
samples produced by our models. Note that there can and should be
disimilarities between these kinematics, which arises from the use of
different gravitational potentials to calculate the kinematic properties, and
varying underlying positional distributions. But overall, on the basis of
kinematic properties, our samples are largely consistent with what is seen in
real observational data.

Another detail we must address is how well these types of DFs represent the
high-$\beta$ stellar halo associated with the debris of GE. This is important
given that we spent a great deal of time in § 5 presenting regions of phase
space that uniquely identified high-$\beta$ samples with high purity, with the
implication that if we applied these to real data they would serve to identify
GE stars. Generalized DFs corresponding to merger debris located in a stellar
halo have not been formulated. This is mostly because the orbit and mass ratio
of the merging galaxy would play a major role in determining the resulting DF.
It is reasonable to assume, however, that the corresponding DF should be
somewhere between that of a thin tidal stream occupying a distinct locus in
action space (e.g. Bovy, 2014), and the spherical DFs we have used here. It
should also be the case that longer timescales and larger mass ratio mergers
tend to make the DF of the debris resemble the latter, rather than the former.
With this in mind we believe it reasonable to assume that the GE merger debris
should resemble a roughly spherical, halo-like DF. This assertion is supported
by a wide range of evidence, including that: the GE debris is distributed over
the whole sky (see Figure 12 in Helmi, 2020), the GE debris can be described
by a triaxial power law density profile (Mackereth & Bovy, 2020), and that
measured stellar ages indicate GE was accreted roughly 10 Gyr ago (Montalbán
et al., 2020), giving the debris ample time to relax.

Additionally, we can look to simulations to further validate our models, and
specifically to support our use of high-$\beta$ halo DFs to model the GE
remnant. Amorisco (2017) presents the kinematic properties of debris from a
number of simulated mergers parameterized by the mass ratio and orbit
circularity of the infalling satellite. They show that the deposited stellar
debris has a power law density profile, with some appearance of a cuttoff
radius around a few halo scale radii, and the degree of radial bias is high.
These facts together support the use of a spherical, radially biased halo DF
to model the remnants. Jean-Baptiste et al. (2017) present a series of 1:10
merger simulations which show that the $E-L_{z}$ distribution of the deposited
stars is broadly consistent with our high-$\beta$ halo models. They show cusps
around $L_{z}=0$ and are roughly bound by extended funnel shaped boundaries.
Most of these mergers have a significant degree of angular momentum to them
though, which manifests in the skewed appearance of the $E-L_{z}$
distribution.

###  6.2 Selection Criteria for Gaia-Enceladus

In § 5, we explored the feasibility of using kinematics alone to separate
high-$\beta$ populations from ergodic and disc contaminants within the
framework of our fiducial model. We found that most kinematic spaces had
regions of high purity, and using well-informed selections, we were able to
achieve simultaneously high purity ($\gtrsim 0.7$) and relatively high
completeness ($\gtrsim 0.2$). Our main goal in providing these assessments of
purity and completeness is to inform and enable identification of GE stars on
the basis of kinematics alone. Chemical and stellar properties of GE
constituent stars are of significant interest currently, as they tell us a
great deal about the progenitor dwarf galaxy. Since it is not obvious a priori
what these properties are in the context of the broader stellar halo, it is
necessary to develop accurate methods to extract GE stars from larger samples
of halo and disc stars. It is for this reason that throughout this paper, we
have placed more emphasis on achieving high purity when it comes to separating
low- and high-$\beta$ populations, than on completeness.

There already exist selection criteria used in the literature to define the
boundaries of GE for chemical and stellar property assays. Here we compare our
results with some of the more common selections. For reference, our benchmark
selection is the union of our action diamond selection ellipse and $e-L_{z}$
selection ellipse (although as we have mentioned this is very similar to the
action diamond selection alone). The most apt comparison for most of the
literature is to our _Survey_ sample, where the purity without (with) disc
contamination is 0.82 (0.81) and completeness of our combined selection is
0.22. In order to compare other selections in a fair manner we evaluate purity
and completeness for them using our fiducial _Survey_ model. We will quote
purity both with and without disc contaminants, but it is important to keep in
mind that disc stars may also be efficiently removed from a sample using
abundances.

First, we consider the selection criterion of Feuillet et al. (2020), which is
also used in Matsuno et al. (2021), and is defined as $\sqrt{J_{R}}>30$ kpc1/2
km1/2 s1/2 and $\lvert L_{z}\rvert<500$ kpc km s-1. The completeness and
purity without (with) disc contaminants of this selection are 0.25 and 0.69
(0.68). Second is the selection of Myeong et al. (2019), also used by Monty et
al. (2020) and Cordoni et al. (2020), which is defined as $\lvert
J_{\phi}/J_{\mathrm{tot}}\rvert<0.07$ and
$(J_{z}-J_{R})/J_{\mathrm{tot}}<-0.3$. The completeness and purity without
(with) disc contaminants of this selection are 0.31 and 0.78 (0.76). A final
selection criterion to consider is that of Naidu et al. (2020), based off of
Belokurov et al. (2018) and similar in spirit to Mackereth & Bovy (2020),
which is defined as $e>0.7$. We note specifically that Naidu et al. (2020)
apply this criteria to their data only after removing a subset of the halo
corresponding to the Sagittarius dwarf and the hot thick disc. The
completeness and purity without (with) disc contaminants of this selection are
0.90 and 0.59 (0.36), as per § 5.2.

Our intent here is not to critique the selection criteria of these studies,
but instead to provide an unbiased characterization of their performance in
the context of our models. What we see from this comparison, and indeed our
study in general, is that when one is kinematically separating high-$\beta$
populations from ergodic populations there is a tradeoff between purity and
completeness. If one’s interest is in maximal purity of a high-$\beta$
population, then the selection criteria we have provided are best,
specifically our action diamond ellipse. By combining this selection with the
$v_{R}-v_{T}$ selection then purity with (without) disc samples is 0.84 (0.83)
and completeness is 0.06. This high-purity selection is almost devoid of
contamination by disc stars, yet has lower completeness. If one is instead
working with a sample with lower numbers of stars and completeness is the
goal, then we recommend our $e-L_{z}$ selection which achieves purity with
(without) disc samples of 0.74 (0.72) and completeness of 0.48.

###  6.3 Structure Identification

One of the most remarkable emergent features of our _Survey_ models was the
appearance of artificial substructure in the distribution of halo samples in
the $E-L_{z}$ plane. The substructure is symmetric about the $L_{z}=0$, and
appears as a small overdensity of stars at low energies. The fact that this
feature appears in our simple halo models defined by a single truncated power
law density profile embedded in a sphericalized Milky Way emphasizes that it
is entirely spurious. We thoroughly explored this phenomenon in § 4.2, and
found it to be caused by a bimodality in the underlying radial distribution of
APOGEE DR16 likely halo stars. Many stars are concentrated near the location
of the Sun, about 8 kpc from the center, and many are concentrated also at
about 3.5 kpc in radius. The cause of this bimodality is twofold. First, many
pointings glance the Galactic bulge (see Figure 1), which means there is a
concentration of halo stars just on its outskirts which share a galactocentric
radius, but have differing galactocentric azimuth. Second, pointings towards
the Galactic bulge contribute the most halo stars right on the outskirts of
the bulge, as opposed to near the Sun, since observed volume scales with the
square of the distance from us, and the density of halo stars should be
centrally concentrated.

This is a particularly pertinent result given that the $E-L_{z}$ plane is
often used as a means of structure identification (e.g. Helmi et al., 2018;
Koppelman et al., 2019b; Horta et al., 2021), and therefore we strongly
recommend that authors consider the selection biases present in their samples
when attempting to verify the nature of newly discovered substructures. This
is especially important for any survey with a non-trivial spatial selection
function, which includes essentially all modern spectroscopic surveys.
Furthermore, while the density fluctuations we observe are symmetric about the
$L_{z}=0$ axis, we have the clarity of being able to separate halo stars from
disc stars. This is not a luxury afforded to most samples, which will also
contain contamination from disc stars. These contaminants would potentially
hide the symmetric nature of structure that arises thanks to the types of
radial selection biases we explore here. This could result in, for example,
the density fluctuations being observed in the lower density retrograde half
of the plane but not the prograde half where they could be masked by disc
contaminants.

We specifically highlight the case of the newly discovered Heracles structure
in the inner Galactic halo by Horta et al. (2021) for two reasons. First, they
also use the APOGEE DR16 data that underpins our models. Second, the structure
they identify shares many features in common with the artifacts we find. This
is to say that Heracles is roughly symmetric about the $L_{z}=0$ axis, and it
appears as a separate overdensity of stars at lower energies than the primary
group of halo stars. We suggest that the structure be further investigated to
confirm that it does not appear as a result of the APOGEE selection function.

In a similar vein, a potentially useful application of the sorts of models we
have presented here, is to use them for the differential identification of
substructures in both the disc and stellar halo. Because our models use
simple, smooth underlying potentials, they can be used to generate phase space
realizations of those smooth potentials which can be compared with
observations. The advantage of this approach is that it would be relatively
simple to account for the types of selection effects we just highlighted,
because phase space need only be sampled at the locations where the real data
are found. A major challenge of this approach is determining the impact of
assuming the nature of the underlying potential, but this could be
circumvented to a degree by comparing many different potentials.

##  7 Conclusion

In this work, we have demonstrated a framework for assessing the kinematic
properties of the types of composite stellar populations which are observed by
modern spectroscopic surveys complemented by Gaia kinematics. Our methodology
relies on a simple, DF-based approach to generate mock kinematics for both
solar-vicinity and APOGEE DR16 data using galpy, and would be simple to
replicate for other surveys. We therefore suggest that these sorts of models
be created for many different surveys in order that their unique kinematic
distributions and underlying selection effects can be properly gauged.

We provide a general overview of many kinematic spaces commonly used in the
study of the Milky Way, and specifically the stellar halo. In each kinematic
space, we show two sets of models: the _Solar_ models which are simple and
mimic stellar samples concentrated near the Sun, and the _Survey_ models which
are more complicated and designed to mimic samples observed with spectroscopic
surveys. Each model uses a stellar halo parameterized by the anisotropy:
$\beta$, which is a useful quantity for describing the two dominant stellar
halo populations in the Milky Way: the low-$\beta$ metal poor halo and the
high-$\beta$ GE remnant.

In each kinematic space we use our models to show where disc populations, low-
and high-$\beta$ halo populations all lie. We use purity as a metric to
specifically assess each space on how well it separates low-$\beta$ halo
populations from high-$\beta$ populations, finding that the scaled action
diamond, the $\sqrt{J_{R}}-L_{z}$ plane, and $e-L_{z}$ plane achieve highest
performance, with a focus on purity. We briefly explored the effects of survey
selection functions on kinematic properties of observed samples, finding that
they can result in the appearance of specious structures.

We summarize this work with a number of recommendations for authors to
consider as they explore current and future datasets:

  * •

The scaled action diamond offers superior performance to other kinematic
spaces when seeking to separate low- and high-$\beta$ stellar halo
populations.

  * •

If invariance to choice of Milky Way potential is most important when
separating low- and high-$\beta$ populations, then consider using the
$v_{R}-v_{T}$ plane, which outperforms the Toomre diagram.

  * •

If high-purity is most important when separating low- and high-$\beta$
populations, then we provide a number of simple-to-implement selection
criteria. We recommend combining our action diamond and $v_{R}-v_{T}$
selections (see Table 2), and we estimate this criterion gives a purity of
0.84 and completeness of 0.06, and is relatively immune to contamination from
disc stars.

  * •

If high completeness is most important when separating low- and high-$\beta$
populations, then we recommend using the $e-L_{z}$ plane. Our selection
achieves purity of 0.74 and completeness of 0.48, and is also relatively
immune to contamination from disc stars.

  * •

In general, the $E-L_{z}$ plane and Toomre diagram do not do as good a job of
separating low- and high-$\beta$ populations. Additionally, eccentricity used
by itself is not a good parameter for selecting these populations. Commonly
used eccentricity cuts for the radially anisotropic halo, such as $e\sim 0.7$
only result in a sample with purity of 0.59. Only extreme eccentricity cuts
($e\gtrsim 0.95$) achieve good purity.

  * •

Be aware of all selection biases in the underlying sample, which can cause
spurious structures to appear in kinematic spaces with explicit dependencies
on radius, specifically in the $E-L_{z}$ plane. We show that in APOGEE DR16, a
bias in the underlying Galactocentric radial distribution of sources results
in the appearance of artifacts in the $E-L_{z}$ plane which could be easily
mistaken for real substructure.

## Acknowledgements

JMML and JB acknowledge financial support from NSERC (funding reference number
RGPIN-2020-04712) and an Ontario Early Researcher Award (ER16-12-061). 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. Funding for the Sloan
Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the
U.S. Department of Energy Office of Science, and the Participating
Institutions. SDSS-IV acknowledges support and resources from the Center for
High-Performance Computing at the University of Utah. The SDSS web site is
www.sdss.org.

## Data Availability

The APOGEE DR16 data used in this article are available at:
https://www.sdss.org/dr16. The Gaia data used in this article are available
at: https://gea.esac.esa.int/archive/.

## References

  * Ahumada et al. (2020) Ahumada R., et al., 2020, ApJS, 249, 3
  * Amorisco (2017) Amorisco N. C., 2017, MNRAS, 464, 2882
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Belokurov et al. (2020) Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020, MNRAS, 494, 3880
  * Bennett & Bovy (2019) Bennett M., Bovy J., 2019, MNRAS, 482, 1417
  * Binney (2010) Binney J., 2010, MNRAS, 401, 2318
  * Binney (2012) Binney J., 2012, MNRAS, 426, 1324
  * Binney & McMillan (2011) Binney J., McMillan P., 2011, MNRAS, 413, 1889
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press, Princeton, NJ
  * Bovy (2014) Bovy J., 2014, ApJ, 795, 95
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Bovy et al. (2012) Bovy J., Rix H.-W., Liu C., Hogg D. W., Beers T. C., Lee Y. S., 2012, ApJ, 753, 148
  * Buder et al. (2020) Buder S., et al., 2020, arXiv e-prints, p. arXiv:2011.02505
  * Conroy et al. (2019) Conroy C., et al., 2019, ApJ, 883, 107
  * Cordoni et al. (2020) Cordoni G., et al., 2020, MNRAS,
  * Cuddeford (1991) Cuddeford P., 1991, MNRAS, 253, 414
  * Cui et al. (2012) Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics, 12, 1197
  * Dalton et al. (2016) Dalton G., et al., 2016, in Evans C. J., Simard L., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. p. 99081G, doi:10.1117/12.2231078
  * Deason et al. (2011) Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903
  * Deason et al. (2019) Deason A. J., Belokurov V., Sanders J. L., 2019, MNRAS, 490, 3426
  * Di Matteo et al. (2019) Di Matteo P., Haywood M., Lehnert M. D., Katz D., Khoperskov S., Snaith O. N., Gómez A., Robichon N., 2019, A&A, 632, A4
  * Erkal et al. (2019) Erkal D., et al., 2019, MNRAS, 487, 2685
  * Fattahi et al. (2019) Fattahi A., et al., 2019, MNRAS, 484, 4471
  * Feuillet et al. (2020) Feuillet D. K., Feltzing S., Sahlholdt C. L., Casagrande L., 2020, MNRAS, 497, 109
  * Forbes (2020) Forbes D. A., 2020, MNRAS, 493, 847
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2020) Gaia Collaboration Brown A. G. A., Vallenari A., Prusti T., de Bruijne J. H. J., Babusiaux C., Biermann M., 2020, arXiv e-prints, p. arXiv:2012.01533
  * Garrow et al. (2020) Garrow T., Webb J. J., Bovy J., 2020, MNRAS, 499, 804
  * Gilmore et al. (2012) Gilmore G., et al., 2012, The Messenger, 147, 25
  * Gravity Collaboration et al. (2019) Gravity Collaboration et al., 2019, A&A, 625, L10
  * Harris (1996) Harris W. E., 1996, AJ, 112, 1487
  * Hawkins et al. (2015) Hawkins K., Jofré P., Masseron T., Gilmore G., 2015, MNRAS, 453, 758
  * Haywood et al. (2018) Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113
  * Helmi (2020) Helmi A., 2020, ARA&A, 58, 205
  * Helmi et al. (1999) Helmi A., White S. D. M., de Zeeuw P. T., Zhao H., 1999, Nature, 402, 53
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Horta et al. (2021) Horta D., et al., 2021, MNRAS, 500, 1385
  * Iorio et al. (2018) Iorio G., Belokurov V., Erkal D., Koposov S. E., Nipoti C., Fraternali F., 2018, MNRAS, 474, 2142
  * Jean-Baptiste et al. (2017) Jean-Baptiste I., Di Matteo P., Haywood M., Gómez A., Montuori M., Combes F., Semelin B., 2017, A&A, 604, A106
  * Kollmeier et al. (2017) Kollmeier J. A., et al., 2017, arXiv e-prints, p. arXiv:1711.03234
  * Koppelman et al. (2018) Koppelman H., Helmi A., Veljanoski J., 2018, ApJ, 860, L11
  * Koppelman et al. (2019a) Koppelman H. H., Helmi A., Massari D., Roelenga S., Bastian U., 2019a, A&A, 625, A5
  * Koppelman et al. (2019b) Koppelman H. H., Helmi A., Massari D., Price-Whelan A. M., Starkenburg T. K., 2019b, A&A, 631, L9
  * Kruijssen et al. (2020) Kruijssen J. M. D., et al., 2020, MNRAS, 498, 2472
  * Lancaster et al. (2019) Lancaster L., Koposov S. E., Belokurov V., Evans N. W., Deason A. J., 2019, MNRAS, 486, 378
  * Leung & Bovy (2019a) Leung H. W., Bovy J., 2019a, MNRAS, 483, 3255
  * Leung & Bovy (2019b) Leung H. W., Bovy J., 2019b, MNRAS, 489, 2079
  * Levi et al. (2019) Levi M., et al., 2019, in Bulletin of the American Astronomical Society. p. 57 (arXiv:1907.10688)
  * Lian et al. (2020) Lian J., et al., 2020, MNRAS, 494, 2561
  * Mackereth & Bovy (2018) Mackereth J. T., Bovy J., 2018, PASP, 130, 114501
  * Mackereth & Bovy (2020) Mackereth J. T., Bovy J., 2020, MNRAS, 492, 3631
  * Mackereth et al. (2019a) Mackereth J. T., et al., 2019a, MNRAS, 482, 3426
  * Mackereth et al. (2019b) Mackereth J. T., et al., 2019b, MNRAS, 489, 176
  * Majewski et al. (2017) Majewski S. R., et al., 2017, AJ, 154, 94
  * Martell et al. (2017) Martell S. L., et al., 2017, MNRAS, 465, 3203
  * Matsuno et al. (2021) Matsuno T., Hirai Y., Tarumi Y., Hotokezaka K., Tanaka M., Helmi A., 2021, arXiv e-prints, p. arXiv:2101.07791
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Miyamoto & Nagai (1975) Miyamoto M., Nagai R., 1975, PASJ, 27, 533
  * Montalbán et al. (2020) Montalbán J., et al., 2020, arXiv e-prints, p. arXiv:2006.01783
  * Monty et al. (2020) Monty S., Venn K. A., Lane J. M. M., Lokhorst D., Yong D., 2020, MNRAS, 497, 1236
  * Myeong et al. (2019) Myeong G. C., Vasiliev E., Iorio G., Evans N. W., Belokurov V., 2019, MNRAS, 488, 1235
  * Naidu et al. (2020) Naidu R. P., Conroy C., Bonaca A., Johnson B. D., Ting Y.-S., Caldwell N., Zaritsky D., Cargile P. A., 2020, ApJ, 901, 48
  * Naidu et al. (2021) Naidu R. P., et al., 2021, arXiv e-prints, p. arXiv:2103.03251
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Necib et al. (2019) Necib L., Lisanti M., Belokurov V., 2019, ApJ, 874, 3
  * Niederste-Ostholt et al. (2010) Niederste-Ostholt M., Belokurov V., Evans N. W., Peñarrubia J., 2010, ApJ, 712, 516
  * Nissen & Schuster (2010) Nissen P. E., Schuster W. J., 2010, A&A, 511, L10
  * Posti et al. (2015) Posti L., Binney J., Nipoti C., Ciotti L., 2015, MNRAS, 447, 3060
  * Sanders (2012) Sanders J., 2012, MNRAS, 426, 128
  * Schönrich et al. (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
  * Steinmetz et al. (2020) Steinmetz M., et al., 2020, AJ, 160, 82
  * Trick et al. (2019) Trick W. H., Coronado J., Rix H.-W., 2019, MNRAS, 484, 3291
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 484, 2832
  * Williams & Evans (2015) Williams A. A., Evans N. W., 2015, MNRAS, 448, 1360
  * Xue et al. (2015) Xue X.-X., Rix H.-W., Ma Z., Morrison H., Bovy J., Sesar B., Janesh W., 2015, ApJ, 809, 144
  * Yanny et al. (2009) Yanny B., et al., 2009, AJ, 137, 4377
  * de Jong et al. (2019) de Jong R. S., et al., 2019, The Messenger, 175, 3
  * van der Marel (2006) van der Marel R. P., 2006, in Livio M., Brown T. M., eds, Vol. 17, The Local Group as an Astrophysical Laboratory. pp 47–71 (arXiv:astro-ph/0404192)