# On the robustness of the velocity anisotropy parameter in probing the
stellar kinematics in Milky Way like galaxies: Take away from TNG50 simulation

Razieh Emami Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Lars Hernquist Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Mark
Vogelsberger Department of Physics, Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Xuejian Shen TAPIR, California Institute of Technology, Pasadena, CA 91125,
USA  Joshua S. Speagle (沈佳士) Banting & Dunlap Fellow  David A. Dunlap
Department of Astronomy & Astrophysics, University of Toronto, 50 St. George
Street, Toronto ON M5S 3H4, Canada  Dunlap Institute for Astronomy and
Astrophysics, University of Toronto, 50 St George Street, Toronto, ON M5S 3H4,
Canada  Department of Statistical Sciences, University of Toronto, 100 St
George St, Toronto, ON M5S 3G3, Canada  Jorge Moreno Department of Physics and
Astronomy, Pomona College, Claremont, CA 91711, USA  Charles Alcock Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  Shy Genel Center for Computational Astrophysics, Flatiron Institute, New
York, USA  Columbia Astrophysics Laboratory, Columbia University, 550 West
120th Street, New York, NY 10027, USA  John C. Forbes Center for Computational
Astrophysics, Flatiron Institute, New York, USA  Federico Marinacci Department
of Physics & Astronomy "Augusto Righi", University of Bologna, via Gobetti
93/2, 40129 Bologna, Italy  Paul Torrey Department of Astronomy, University of
Florida, 211 Bryant Space Sciences Center, Gainesville, FL 32611, USA

(Accepted in ApJ)

###### Abstract

We analyze the velocity anisotropy of stars in real and energy space for a
sample of Milky Way-like galaxies in the TNG50 simulation. We employ different
selection criteria, including spatial, kinematic and metallicity cuts, and
make three halo classes ($\mathcal{A}$-$\mathcal{C}$) which show mild-to-
strong sensitivity to different selections. The above classes cover 48%, 16%
and 36% of halos, respectively. We analyze the $\beta$ radial profiles and
divide them into either monotonically increasing radial profiles or ones with
peaks and troughs. We demonstrate that halos with monotonically increasing
$\beta$ profiles are mostly from class $\mathcal{A}$, whilst those with
peaks/troughs are part of classes $\mathcal{B}$-$\mathcal{C}$. This means that
care must be taken as the observationally reported peaks/troughs might be a
consequence of different selection criteria. We infer the anisotropy parameter
$\beta$ energy space and compare that against the $\beta$ radial profile. It
is seen than 65% of halos with very mild sensitivity to different selections
in real space, are those for which the $\beta$ radial and energy profiles are
closely related. Consequently, we propose that comparing the $\beta$ radial
and energy profiles might be a novel way to examine the sensitivity to
different selection criteria and thus examining the robustness of the
anisotropy parameter in tracing stellar kinematics. We compare simulated
$\beta$ radial profiles against various observations and demonstrate that, in
most cases, the model diversity is comparable with the error bars from
different observations, meaning that the TNG50 models are in good overall
agreement with observations.

Milky Way Galaxy, TNG simulation, star, velocity anisotropy profile,
eccentricity, prograde/radial/retrograde orbits, metallicity

##  1 Introduction

Halo stars comprise only $\sim 1\%$ of the stellar mass in the Milky Way (MW).
Nevertheless, they provide remarkable information about the assembly history
of our Galaxy. Stars belonging to the stellar halo are typically old and
metal-poor, pointing to the halo’s ancient origin. By tracking the current
orbits of stars, we may get some information on their motions at the early
times. Consequently, we may use their current orbit to track the kinematics of
their progenitors, such as external satellite galaxies or the gas clouds where
these stars formed in initially (Hattori et al., 2013). Furthermore, stellar
kinematics can also be employed for studying the distribution of the gaseous
component (Cinzano & van der Marel, 1994; Krajnović et al., 2005; Adams et
al., 2012, 2014).

A widely used technique to model the galactic mass distribution is the Jeans
dynamical modelling (Jeans, 1915; Binney, 1980; Merritt, 1985; Dejonghe &
Merritt, 1992; Courteau et al., 2014). Motivated by this, it is customary to
compute the root-mean-square (rms) of the radial velocity as $v_{\rm
rms}\equiv\sqrt{\langle v^{2}_{r}\rangle}$ (El-Badry et al., 2017) and the
velocity anisotropy parameter,

$\beta\equiv
1-\left(\frac{\sigma^{2}_{\theta}+\sigma^{2}_{\phi}}{2\sigma^{2}_{r}}\right).$
(1)

Here $\sigma_{i}$ refers to the velocity dispersion along the following
directions: $i=(r,\theta,\phi)$, and it is calculated as
$\sigma_{i}\equiv\sqrt{\langle v^{2}_{i}\rangle-\langle v_{i}\rangle^{2}}$.
The velocity anisotropy parameter, hereafter $\beta$, was initially introduced
by Binney (1980) to measure the orbital structure of a given system. In this
parametrization, $\beta$ ranges from $-\infty$ (for purely tangential orbits)
through $0$ (for completely isotropic motions) and to $1$ (for purely radial
orbits). Radial orbits correspond to $\beta>0$, whilst tangential ones have
$\beta<0$. The velocity anisotropy is vital in accurately measuring the mass
distribution in dispersion-supported systems, from dwarf spheroidal galaxies
(Gilmore et al., 2007; Strigari et al., 2007; Wolf et al., 2010) to massive
ellipticals (Dekel et al., 2005).

Direct measurements of $\beta$ have been somewhat challenging because they
require a full three-dimensional map of stellar velocities, which had not been
possible until very recently, thanks to the development of Hubble Space
Telescope (HST) and the Gaia mission. Observers commonly measure the line-of-
sight (LOS) velocity and use this to infer the $\beta$ (see for example Sirko
et al., 2004; Kafle et al., 2012; Deason et al., 2012; King et al., 2015, and
references therein). The main reason for this is that owing to our position in
the MW galaxy, line-of-sight (LOS) velocities provide useful information about
tangential velocities. However, Hattori et al. (2017) pointed out that since
stellar orbits become progressively more radial at larger galactocentric
distances, using only LOS velocities may underestimate the $\beta$.
Consequently, the velocity anisotropy parameter has instead only been measured
indirectly via dynamical modelling (see for example Cinzano & van der Marel,
1994; Diakogiannis et al., 2014a, b, and references therein).

Fortunately, it is now possible to measure tangential velocities for a small
sample of stars away from the galactic center by using full three-dimensional
velocity measurements (e.g. Cunningham et al., 2016; Smith et al., 2009; Bond
et al., 2010; Sohn et al., 2012; Watkins et al., 2019; Gaia Collaboration et
al., 2018, and references therein), although getting very accurate LOS
velocity may require obtaining the spectra.

In observational studies, authors have adopted different assumptions for the
$\beta(r)$ profile. For example, it is commonly assumed that either
$\beta(r)$= Const (e.g., Łokas et al., 2005; Koch et al., 2007; Walker et al.,
2007; Battaglia et al., 2008; Łokas, 2009; Walker et al., 2009, and references
theirin), or that $\beta(r)$ varies radially (e.g. Kleyna et al., 2001;
Wilkinson et al., 2004; Mamon & Łokas, 2005; Gilmore et al., 2007; Battaglia
et al., 2008; Mamon et al., 2013; Mashchenko, 2015). Whilst a constant $\beta$
is well suited to measure the mass of the MW, a radially-varying $\beta$, may
have additional power in probing the MW’s accretion history (e.g. Cunningham
et al., 2019; Loebman et al., 2018; Bird et al., 2020). Morrison et al. (1990)
used the G and K giants (located up to few kpc from the Sun) and estimated
$\beta\approx 0.5$. Chiba & Yoshii (1998) estimate $\beta=0.52\pm 0.07$ from a
sample of metal-poor halo stars near the Sun. Beyond the Solar neighborhood,
one may also estimate $\beta$ with line-of-sight velocities, by adopting mass-
distribution models of the galaxy. Appropriate tracers include blue horizontal
branch (BHB) stars, K giants and F-type stars (see e.g. Sommer-Larsen et al.,
1994; Thom et al., 2005; Kafle et al., 2012; King et al., 2015, and references
therein). Smith et al. (2009) used halo sub-dwarfs from the Sloan Digital Sky
Survey (SDSS) and estimated $\beta=0.69\pm 0.01$. Bird et al. (2019) used K
giants from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope
(LAMOST) catalog and measure $\beta$ as a function of Galactocentric radius.

On the theory front, simulations of MW-like galaxies often result in radially-
increasing $\beta$ profiles (e.g. Diemand et al., 2005; Abadi et al., 2006;
Sales et al., 2007; Rashkov et al., 2013; Kafle et al., 2012; Stinson et al.,
2013). As shown by Loebman et al. (2018), analyses of this kind unveil more
complex structures in $\beta$ and in the halo in general (e.g. Beers et al.,
2012; Schuster et al., 2012; Zuo et al., 2017).

In this paper, we analyze the stellar kinematics inferred from the stellar
velocity anisotropy parameter, $\beta$, in a sample of Milky Way (MW) like
galaxies from the TNG50 run of the Illustris-The-Next-Generation (Illustris-
TNG) project. The enhanced resolution in TNG50 makes it possible to analyze
radial $\beta$ profiles very accurately, near the level of the zoom-in
simulations and also consider the impact of cosmic evolution and mergers of a
larger statistical population as TNG50 simulation has many more galaxies in
the box. We analyze the $\beta$ radial and energy profiles at redshift $z=0$
and investigate their dependencies on various selection criteria, including
spatial and kinematic cuts as well as the ones coming from some metallicity
cuts.

We design a halo classification based on the sensitivity levels of $\beta$
profile to various selection criteria and place halos in three classes.
$\mathcal{A}$ class refers to halos with very mild sensitivity to different
selections. Halos part of class $\mathcal{B}$ show some levels of
sensitivities to different selections and those associated with $\mathcal{C}$
show a very strong dependencies to model selections. We show that 48%, 16% and
36% of halos are part of classes $\mathcal{A}$ to $\mathcal{C}$, respectively.
We analyze the $\beta$ radial profile and show that it can be divided to
monotonically increasing profiles or the ones with peaks/troughs, which
contain 32% and 68% of the halos, respectively. It is shown that, the $\beta$
profile from the monotonically increasing class smoothly increases from the
interior to the exterior part of the halo. While, the $\beta$ from the class
with peaks/troughs experiences some local fluctuations. We further demonstrate
that almost all of the halos from the first class are part of category
$\mathcal{A}$, with only one exception. On the contrary, members of the class
with peaks/troughs are mostly part of category $\mathcal{B}$ and
$\mathcal{C}$, with only 29% of them being part of category $\mathcal{A}$.
This means that care should be taken in interpreting the observed
peaks/troughs as they might depend on the actual selection criteria.

We infer the impact of different metallicity and eccentricity cuts in the
$\beta$ radial profile and also investigate the impact of different stellar
orbital types such as the radial, prograde and retrograde on the $\beta$
radial profile.

Moving in the energy space, we infer and compare the $\beta$ radial and energy
profiles against each other and make another halo classification based on
their differences. It is seen that 65% of the halos with very mild sensitivity
on different selection criteria, in the real space, are among those for which
the $\beta$ radial and energy profiles are closely related. On the contrary,
halos in which the $\beta$ radial and energy profile are rather different are
entirely coming from category $\mathcal{C}$. This again establishes that one
should be careful when drawing any conclusions from these classes.
Furthermore, we propose that one quantitative way to check if the
observational results are very robust is to also compute the $\beta$ energy
profile and compare that with the $\beta$ radial profile. Based on our
conclusion it is likely that those with very similar profiles are also very
robust against altering the selection criteria.

It must be emphasized that our current $\beta$ analysis in the real and energy
spaces, and their correlations, are purely theoretical and are based on a full
awareness of the gravitational potential for the MW. From the observational
perspective, we may only estimate the potential to some extent. Furthermore,
the current paper is not aimed to prove that we identify an unbiased sample of
stars using this comparison. But rather that for a given sample of stars such
a comparison between the $\beta$ radial and energy profiles likely indicates
how the sample may show sensitivities to changing different selection
criteria.

Motivated by the aforementioned halo classification, throughout this paper, we
make a sub-halo sample of 4 galaxies from each of the above classes and study
their behavior thoroughly. This includes halos [4, 12, 20, 25] which are the
members of Category [$\mathcal{A}$, $\mathcal{B}$, $\mathcal{B}$,
$\mathcal{C}$], respectively.

Ultimately, we overlay our theoretical outcomes, with few different selection
criteria, on top of the most recent observational results and compare them
with each other. Our comparison is made both at the level of isolated data
points as well as the extended observations in which the $\beta$ profile is
reported in more than one single point. It is shown that there are good levels
of agreement between the two both for the isolated points and for the extended
observations where we see for each halo there are some data points that easily
pass through the profile.

This paper is organized as follows. Section 2 summarizes the TNG50 simulation
and our sample of MW-like galaxies. Section 3 links stellar metallicity to
other stellar properties. Section 4 computes the radial and energy velocity
anisotropy profiles. Section 5 compares our simulated radial-$\beta$ profiles
against observations. Lastly, Section 6 presents our conclusions.

##  2 TNG50 simulation

TNG50 is the highest resolution run of the suite of large-scale series of
IllustrisTNG simulations (e.g. Pillepich et al., 2019; Nelson et al., 2019a).
It provides an excellent combination of the volume and resolution close to the
level of zoom-in simulations (e.g., Monachesi et al., 2016; Grand et al.,
2018; Hani et al., 2019) and thus provides a natural avenue to investigate the
impact of galaxy evolution with very high resolution. It evolves supermassive
black holes (SMBHs), dark matter, gas, stars and magnetic fields within a
periodic boundary volume of 51.7 kpc3. Its softening length is 0.39 comoving
kpc $h^{-1}$ for $z\geq 1$, reducing to 0.195 proper kpc $h^{-1}$ for $z<1$.

In TNG50, the model for galaxy formation includes a stochastic, gas-density-
threshold-based for the star formation, the evolution of mono-age stellar
populations which are being represented by the star particles, the chemical
enrichment of the Inter Stellar Medium (ISM) as well as the tracking of 9
different chemical elements including (H, He, C, N, O, Ne, Mg, Si, Fe) in
addition to total gas metallicity and the Europium. It also includes the gas
cooling and heating, the feedback from supernovae which is in the form of
galactic winds, seeding and the growth of supermassive black holes (SMBHs) and
the injection of the energy as well as the momentum from SMBHs into the
surrounding gas.

Using various selection criteria, such as dark matter halo mass and choosing
the rotationally-supported stellar population in (Emami et al., 2020, 2021a,
2021b), we created a sample of 25 MW-like galaxies in the TNG50 simulation.
This paper probes the structure of the stellar distribution, ignoring the
impact of satellites, with stellar kinematics as the tracer. We also analyze
different stellar properties, including stellar age, metallicty and the
stellar velocity anisotropy profile. One of our goals is to investigate how
different selection criteria based on these quantities impact radial velocity
anisotropy profiles.

Table 1: Linking between the subhalo ID with the galaxy number in our sample of MW like galaxies in TNG50. $1\mapsto 476266$ | $2\mapsto 478216$ | $3\mapsto 479938$ | $4\mapsto 480802$ | $5\mapsto 485056$
---|---|---|---|---
$6\mapsto 488530$ | $7\mapsto 494709$ | $8\mapsto 497557$ | $9\mapsto 501208$ | $10\mapsto 501725$
$11\mapsto 502995$ | $12\mapsto 503437$ | $13\mapsto 505586$ | $14\mapsto 506720$ | $15\mapsto 509091$
$16\mapsto 510585$ | $17\mapsto 511303$ | $18\mapsto 513845$ | $19\mapsto 519311$ | $20\mapsto 522983$
$21\mapsto 523889$ | $22\mapsto 529365$ | $23\mapsto 530330$ | $24\mapsto 535410$ | $25\mapsto 538905$
Figure 1: The probability density distribution of stellar metallicity in our
MW-like galaxy sample. It is evident that stars from different galaxies have
quite similar metallicity distribution. Different colors refer to various
galaxies. Low metallicities are in the tail of the distribution.

##  3 Stellar Metallicity

Since the metallicity cut is one of the key ingredients in our stellar sample
selection, in what follows we study the correlation between stellar
metallicity and various other stellar properties, including the stellar age,
the spatial metallicity distribution and the stellar velocity.

First, we present the stellar metallicity distribution. Figure 1 presents the
1D probability density distribution of all of stars in different galaxies in
our sample. Quite interestingly, stars from different galaxies in our sample
seem to have very similar metal distribution. Moreover, it is evident that
stars with low metallicty are in the tail of the distribution, whilst those
with [Fe/H] in the range [-0.5-0.0], sit on the peak of the distribution. It
is intriguing to see how different metallicty-based cuts might affect the
velocity anisotropy of stars. In Sec. 4.2, we get back to this question in
depth.

Figure 2: 2D Correlation of stellar age vs metallicity for a sub-sample of
four galaxies from our MW-like galaxy sample. The Spearman correlation is
negative indicating that metal-poor stars are older than the metal-rich ones.

###  3.1 Stellar age vs metallicity

Next, we study the correlation between the stellar age and the stellar
metallicity. Figure 2 presents $\mathrm{[Fe/H]}$ versus stellar age for a sub-
sample of four galaxies from our sample, including galaxies [4, 12, 20, 25].
As already stated above, these halos are members of distinct categories. In
each case, we have inferred the Spearman correlation, hereafter r, between the
stellar metallicity and the age.

From the plot, it is inferred that while in some cases, such as galaxy 4, the
stellar metallicity is smoothly declining with increasing stellar age, others
such as galaxy 12 show a break and thus a discrete profile of age vs [Fe/H]
indicating that their associated stars might have different origins. For
instance, they might be accreted as a satellite merged onto the galaxy or have
arisen from star forming gas out of a galaxy merger etc. Intriguingly, such a
break reduces the Spearman correlation as we end up having a bi-modal
age-[Fe/H] relation.

Furthermore, it is seen that in some cases, such as galaxy 20, there is a peak
of star formation at an intermediate time, around 4.5 Gyr, that might be an
indication of a major merger in the history of the galaxy. Comparing the
Spearman correlation of galaxy 4 and 20, we may argue that such a local peak
may have slightly diminished the r in galaxy 20. This is better seen when we
compare these with galaxy 25, in which we have few distinct bursts of stars
and that its |r| is slightly higher than galaxy 20 while it is still less than
the one for galaxy 4.

In summary, our analysis demonstrate that in all cases, $|r|\geq 0.5$, showing
a moderate to large negative correlation between the stellar age vs
metallicity depending on galaxy ancient history. Owing to the above
correlation between the stellar age vs metallicty, different metallicity cuts
naturally pick stars with different ages and allow us to indirectly probe the
velocity distribution of young vs old stars. In Sec. 4.2, we address this
question.

Figure 3: 2D Correlation of the stellar distance vs stellar metallicity for a
sub-sample of four galaxies from our MW-like galaxy sample. The Spearman
correlation is modest and negative indicating metal-rich stars are located
closer to the galactic center.

###  3.2 Spatial distribution of metallicity

Below we explore the possible correlation between the location of stars and
the stellar metallicity. Figure 3 presents the 2D distribution of the stellar
radii vs the stellar metallicity for a sub-sample of 4 out of our galaxy
sample. From the plot, it is inferred that metal-rich stars are generally
closer to the galaxy center than the metal-poor ones, in line with the
observations (e.g. Huang et al., 2015; Anders et al., 2017; Santucci et al.,
2020; Dietz et al., 2020). Combining that with the above age-metallicity
relation, we argue that galaxies in our sample mostly have inside-out star
formation in agreement with the recent studies of the inside-out growth of
stellar disk (see e.g. Frankel et al., 2019; Bird et al., 2021; Johnson et
al., 2021, and references therein). The slope of the metallicity-gradient
varies from galaxy to galaxy. Some galaxies show a more extended metallicity
gradient whilst others have a more concentrated profiles. Furthermore, the
extent to which we infer the dominant population of stars is also correlated
with their age-metallicity profile. For instance, galaxy 12 shows a more
extended profile inline with our former argument of being affected by mergers.
This makes sense as galaxy mergers are expected to expand and dis-locate stars
from their original place. Finally, the Spearman correlation between the
stellar location vs the stellar metallicity is relatively low, indicating that
there are some other factors in locating the stars than their metallicity.

Figure 4: 2D Correlation of the stellar velocity amplitude vs the stellar
metallicity for a sub-sample of four galaxies from our MW-like galaxy sample.
It is inferred that the stellar velocity is largely uncorrelated with the
stellar metallicity as stellar speed is mainly concentrated between 200-300 km
sec-1 in agreement with the observed rotational velocity of stars in our Milky
Way galaxy.

###  3.3 Metallicity vs velocity

As a final step, here we study the correlation between the stellar speed vs
the stellar metallicity. Figure 4 presents the 2D distribution of the stellar
velocity amplitude vs the stellar metallicity for a sub-sample 4 out of our
galaxy sample. In each case, we infer the Spearman correlation between these
quantities demonstrating that they are uncorrelated. This makes sense as our
galaxy sample obeys the well-known Tully-Fisher relation (Tully & Fisher,
1977), as well as the mass-metallicity relation (Gallazzi et al., 2005; Kewley
& Ellison, 2008). The velocity peaks stand between 200-300 km sec-1,
indicating that the kinematic of stars inferred in this sample of MW like
galaxies are consistent with the rotational velocity of stars seen in Milky
Way disk.

Figure 5: $\beta$ radial profile for stars inferred using different selection
criteria. Dashed-black, blue-dotted-dashed, solid-red and dotted-yellow lines
present $\beta$ inferred from stars with kinematic, metallicity based, spatial
cut and no-cut as was used in computing the $\beta$, respectively.

##  4 Velocity anisotropy: $\beta$

In this section, we analyse the radial profile of the velocity anisotropy as
defined in Eq. (1). In our calculations, we use various selection criteria.
Recall that the $\beta$ parameter characterizes various stellar orbital types.
We present the 2D distribution map for the stellar location out of different
selection criteria. Furthermore, we study the impact of changing the
metallicity cuts, eccentricity and different orbital types on $\beta$ radial
profile, $\beta(r)$. Finally, we analyse the $\beta$ energy profile and
compare that with the $\beta$ radial profile.

###  4.1 The $\beta$ radial profile

To infer the $\beta(r)$ for each galaxy, we adopt a coordinate system where
the z-axis (oriented along with the total angular momentum of stars) is
orthogonal to the disk, with (x-y) being chosen randomly (right-handed and
orthogonal to each other) on the disk plane. As a next step, since $\beta(r)$
deals with radial and tangential vectors, we make a transformation to the
spherical coordinate system (for both of the coordinate and the velocity
vectors). We adopt 15 linearly-distributed radial bins between 1 and 100 kpc
and, for each of them, we compute the radial and tangential velocity
dispersions. We force each bin to contain at least 100 stellar particles (in
practice, many of them have over 1000 stellar particles). We choose different
selection criteria for stars. Figure 5 presents the $\beta$ radial profile for
4 different cases as listed below:

$\bullet$ The dashed-black-line refers to the kinematically-chosen stars with
$\epsilon<0.7$, where $\epsilon$ is defined by Eq. (3) in Emami et al.
(2021a). Theoretically, this means that we are choosing halo stars (see e.g.
Monachesi et al., 2019, and references therein). We emphasize that the above
threshold is not unique. Nevertheless, it is robust to small variations.

$\bullet$ The blue-dotted-dashed-line describes a hybrid selection of
$\mathrm{[Fe/H]}<-1.3$ with $|Z|\geq 5$ kpc. This choice is based on the
LAMOST data (Bird et al., 2019).

$\bullet$ With the solid-red-lines, we relax the metallicity cut and only
choose stars above the disk plane $|Z|\geq 5$ kpc.

$\bullet$ The dotted-yellow-line refers to the case with no spatial cut.

At the first glance, it seems that in some cases the velocity anisotropy of
stars with completely different selection criteria might be somewhat close.

Before we analyse the $\beta$ radial profile above in depth, here we present a
spatial map of stars inferred from different selection cuts. This enables us
to develop an intuition on the impact of different selection criteria on the
spatial distribution of stars. Figure 6 presents an image for a typical
galaxy, galaxy 1, along with three different projections. Top to bottom rows
present the 2D surface density map of stars with $\tilde{r}\geq 1$,
$\tilde{r}\geq 1,\epsilon\leq 0.7$, $\tilde{r}\geq 1,\tilde{Z}\geq 5$ and
$\tilde{r}\geq 1,\tilde{Z}\geq 5,$ [Fe/H] $\leq$ -1.3, respectively. Here
$\tilde{X}\equiv X/$kpc with $X=(r,Z)$. Interestingly, the kinematic cut of
$\epsilon\leq 0.7$ and the spatial cut of $\tilde{Z}\geq 5$ both decrease the
spiral structure of the galaxy in the disk plane. While, the current
metallicity cut diminishes the population of stars substantially, as [Fe/H]
$\leq-1.3$ sits on the tail of the stellar distribution.

Figure 6: Face-on (left) and (two) edge-on (middle-right) projections of stars
in galaxy 1 in our sample. From top to bottom rows in the image present the
spatial distribution of stars from $\tilde{r}\geq 1$, $\tilde{r}\geq
1,\epsilon\leq 0.7$, $\tilde{r}\geq 1,\tilde{Z}\geq 5$ and $\tilde{r}\geq
1,\tilde{Z}\geq 5,$ [Fe/H] $\leq$ -1.3, respectively. The color bars refer to
the surface number density of stars in different projections.

Having fully described the details of different lines as presented in Figure
5, below we analyze the $\beta$ radial profile inferred from each of them in
depth and make a halo classification based on the behavior of $\beta$ radial
profile.

####  4.1.1 $\beta$ inferred halo classification from selection criteria

Having presented the $\beta$ radial profile for halos in our sample, here we
make a halo classification based on the sensitivity of $\beta$ radial profiles
to different selection criteria. The goal is to quantify the sensitivity of
$\beta$ to different stellar orbital types achieved from different selections.

In our analysis, for every galaxy, we use the following empirical classifier:

$\delta\beta\equiv\frac{\sum_{i}|\beta_{i}^{\mathrm{Fe}}-\beta_{i}^{\mathrm{No-
cut}}|}{\sum_{i}\mid\beta_{i}^{\mathrm{No-cut}}\mid}.$ (2)

where $\beta_{i}\equiv\beta(r_{i})$ is used for the sake of brevity. In
addition, in Eq. (2), $\beta^{\mathrm{Fe}}_{i}$ refers to the hybrid selection
of $\mathrm{[Fe/H]}<-1.3$ with $|Z|\geq 5$ kpc, while $\beta_{i}^{\mathrm{No-
cut}}$ stands for the case with no spatial cut, as described above. The
selection of $\beta^{\mathrm{Fe}}_{i}$ and $\beta_{i}^{\mathrm{No-cut}}$ in
the halo classifier is owing to the visual differences seen in the curves
associated with them, as presented in Figure 5.

Figure 7: The radial profile of
$\Delta\beta(r)\equiv\left(\beta^{\mathrm{Fe}}/\beta^{\mathrm{No-
cut}}\right)(r)-1$ for halos from different categories. As stated in Section
4.1.1, halos from Categories $\mathcal{A}$ through $\mathcal{C}$ show little
to very strong response on different selection criteria for $\beta$,
respectively. In addition, it is inferred that on average, halos from category
$\mathcal{A}$ show smoother radial profile than the other two classes. This is
more clarified in Figure 8.

As a next step, we define distinct halo classes based on the magnitude of
$\delta\beta$ as listed below:

$\bullet$ Category $\mathcal{A}$: $\delta\beta\leq 0.35$. Halos part of this
category show the least sensitivity to different selection criteria. There are
12 halos in this category, including halos [1,4,5,6,7,13,15,16,17,18,22,23].
This shows $\beta$ is a very robust quantity in about 48% of halos in our
sample.

$\bullet$ Category $\mathcal{B}$: $0.35\leq\delta\beta\leq 0.70$. Halos in
this class, show a mild response to different selection criteria. There are 4
halos in this class, including halos [10, 12, 21, 25].

$\bullet$ Category $\mathcal{C}$: $\delta\beta\geq 0.70$. Halos belong to this
class show a strong dependency to different selections in $\beta$ calculation.
There are in total 9 halos in this class, including halos
[2,3,8,9,11,14,19,20,24]. This demonstrates that in 36% of halos in our
sample, $\beta$ is not a robust quantity as it depends on different selection
criteria.

Figure 7 presents the radial profile of the fractional difference between
$\beta^{\mathrm{Fe}}_{i}$ and $\beta_{i}^{\mathrm{No-cut}}$ for halos from
different classes. In each case, we present:

$\Delta\beta(r)\equiv\left(\frac{\beta^{\mathrm{Fe}}(r)}{\beta^{\mathrm{No-
cut}}(r)}\right)-1.$ (3)

eliminating the points for which the denominator gets below 0.05. From the
plot, it is inferred that halos belong to category $\mathcal{A}$ in average
show smoother radial profile than those from $\mathcal{C}$. Below, we make
this point more clear.

####  4.1.2 Detailed radial profile of $\beta$

Having identified various halo classes, in what follows, we study the $\beta$
radial profile in some depth and make a connection to the aforementioned halo
categories. Considering different profiles, we are dealing with the following
behaviors:

Figure 8:  (Left panel) the profile of $\beta$ with monotonically increasing
radial profile. 32% of halos are in this category. Almost all of them are from
category $\mathcal{A}$, with one exception belong to class $\mathcal{B}$.
(Right panel) the $\beta$ radial profile for cases with peaks and troughs.
About 68% of halos are belonged to this category from which only 29% of halos
are belong to category $\mathcal{A}$. This means that care must be taken in
interpreting the observed peaks/troughs in the $\beta$ radial profile as it
might be an indication of different selection criteria.

$\bullet$ Monotonically increasing: From Figure 5, it is inferred that in 32%
of halos in our sample, $\beta$ increases with the radius implying that
farther out from the center, the orbits become progressively more radial
(i.e., less tangential) and thus less rotationally-dominated. This is in
agreement with e.g. El-Badry et al. (2017). In such cases, $\beta(r)$ profiles
are monotonically increasing (i.e., galaxies 6, 7, 13, 15, 18, 22, 23 and 25).
There are no visible peaks or troughs in the $\beta$ radial profile. Such
galaxies are also not sensitive to our various selection criteria (the various
lines overlap). So $\beta$ radial profile is less biased at identifying
distinct stellar orbital types. Comparing the halos from this class with the
above halo categories, it is evident that almost all of the halos in this case
are from category $\mathcal{A}$. The only exception to this is halo 25 which
is part of class $\mathcal{B}$.

$\bullet$ Profile with Peaks/Troughs: On the contrary, in 68% of cases (e.g.,
galaxies 1, 2, 3, 4, 5, 8, 9, 10, 11, 12, 14, 16, 17, 19, 20, 21 and 24),
$\beta(r)$ experience visible peaks and troughs, in line with the recent
observations (e.g. Kafle et al., 2012; Cunningham et al., 2019). About 29% of
halos in this class are belong to category $\mathcal{A}$ while the rest are
from classes $\mathcal{B}$ and $\mathcal{C}$. As the vast majority of halos
with peaks/troughs show mild/strong response to different selection criteria,
care is warranted because the presence or absence of such peaks and troughs
might actually depend on different selection criteria. For instance, selecting
stars with $\mathrm{[Fe/H]}<-1.3$ typically results in $\beta(r)$ profiles
without these features. Figure 8 compares the details of $\beta$ profiles in
each of the aforementioned cases. On the left panel, we show the $\beta$
radial profile for monotonically increasing cases, while on the right panel,
we present the cases with peaks/troughs.

###  4.2 Impact of different metallicity cuts on $\beta$ radial profile

Having presented the $\beta$ radial profile with the LAMOST based metallicity
cut, here we generalize such a selection criteria and study the impact of
changing different metallicity cuts on the $\beta$ radial profile. As already
stated in Sec. 3.1, since the stellar metallicity is negatively correlated
with the stellar age, various metallicity cuts somewhat choose stars with
different ages. It is then very intriguing to study the impact of the stellar
age on the velocity anisotropy.

Figure 9 presents the $\beta$ radial profile for stars with different spatial
cuts vs various metallicity criteria. The upper panel presents the $\beta$
profile for stars located at $r\geq 1$ kpc, while the bottom one shows the
$\beta$ profile for stars located at $r\geq 1$ kpc and $Z\geq 5$ kpc. In each
row, from the left to the right, we expand over the metallicity cut, from
[Fe/H] $\leq-1.3$, -1.3 $\leq$ [Fe/H] $\leq$ 0 to [Fe/H] $\leq$ 0.0,
respectively. In line with the conclusion of Sec. 3.1, we might expect that
chosen stars in the left panel are on average older than those in the middle
one, while stars on the right panel are the mixture of the young and old
stars. It is inferred that the left panels have higher $\beta$s than the
middle ones. This makes sense, as old stars are less rotationally supported
implying lower $\sigma_{\theta}$ and $\sigma_{\phi}$ and thus higher $\beta$.
On the contrary, young stars have higher tangential velocity components which
correspond to lower values of the $\beta$. Finally, mixed stars, on the right
panel, have $\beta$s in between the young and old populations. Interestingly,
such a conclusion is almost independent of the spatial cut, as the top and
bottom rows are generally quite similar.

Figure 9:  $\beta$ radial profile for stars located at $r\geq 1$ kpc (first
row) as well as those with $r\geq 1$ kpc $\&Z\geq 5$ kpc (second row). In each
row, from the left to the right, we consider stars with [Fe/H] $\leq-1.3$,
-1.3 $\leq$ [Fe/H] $\leq$ 0 and [Fe/H] $\leq$ 0.0. While in the left panels we
present the $\beta$ itself, to facilitate the comparison, in the middle and
right panels, we compute the $\delta\beta$ subtracting the $\beta$ from the
left panel with [Fe/H] $\leq-1.3$. It is clearly seen that expanding over the
range of metallicity, on-average, diminishes the $\beta$.  Figure 10: $\beta$
profile for stellar on different orbital types. From the left to right, we
study the $\beta$ radial profile for all, retrograde, radial and prograde
orbits, respectively. It is generally seen that stars on radial orbits have
slightly higher $\beta$ than the retrograde and the prograde orbits. On the
contrary, stars on retrograde and prograde orbits show similar trends in
$\beta$ radial profile.

###  4.3 Impact of different orbital types on the $\beta$ radial profile

Next, we explore the impact of stellar orbital types on the $\beta$ radial
profile. Figure 10 shows $\beta(r)$ profiles for the retrograde, radial and
prograde orbits (plus the entire stellar set on the left-hand panel). In each
case, we compute $L^{i}_{z}\equiv\vec{L}^{i}\cdot\vec{\hat{L}}_{\mathrm{tot}}$
for $i$-th star, where $\vec{L}^{i}$ refers to the angular momentum of the
$i$-th star while the $\vec{\hat{L}}_{\mathrm{tot}}$ describes the unit vector
along with the total angular momentum of stars. We define our three stellar
orbital types as follows (Naidu et al., 2020):

$a.$ Retrograde orbits: $L^{i}_{z}\leq-500~{}\mathrm{kpc\times
km~{}sec^{-1}}$,

$b.$ Radial orbits: $-500\leq L^{i}_{z}/\left(\mathrm{kpc\times
km~{}sec^{-1}}\right)\leq 500$,

$c.$ Prograde orbits: $L^{i}_{z}\geq 500~{}\mathrm{kpc\times km~{}sec^{-1}}$.

While there are some variations on the exact behavior of individual galaxies,
it is generally true that stars on radial orbits have slightly higher $\beta$
than the retrograde and the prograde orbits. That is expected as radial orbits
have generally higher $\sigma_{r}$ and thus their $\beta$ would be higher
based on Eq. (1). On the contrary, stars on retrograde and prograde orbits
seem to have similar trends though stars on prograde orbits have slightly
higher $\beta$s than the retrograde ones.

Figure 11: $\beta$ profile for different eccentric orbits. From the left to
right we present the entire of stars, stars with $e\leq 0.3$, $0.3\leq e\leq
0.7$ and $e\geq 0.7$, respectively. It is noted that low-eccentricity stars
have lower $\beta$ than the parent population as such stars are on nearly-
circular orbits, On the other contrary, stars with higher eccentricities are
more on radial orbits than tangential, producing higher $\beta$ values.
Figure 12: 1D distribution of the eccentricity of stars in our galaxy sample.
Quite interestingly, most galaxies have an eccentricity peak between 0.2-0.4.
In very few cases, there is a second peak at higher eccentricities of about
1.0 coming from those stars located at larger distances.

###  4.4 Impact of eccentricity on the $\beta$ radial profile

As the next case, here we analyze the impact of the orbital eccentricity on
the $\beta$ radial profile. Where the eccentricity is approximately defined
as,

$e\simeq\sqrt{1-\left(\frac{L}{L_{\mathrm{circ}}}\right)^{2}},$ (4)

where $L_{\mathrm{circ}}\simeq\sqrt{GM(r_{c})r_{c}}$ with $M(r_{c})$ referring
to the total mass interior to $r_{c}$ describing the circular radii, see
Equation 3 in Merritt (2015) and Emami et al. (2021a) for more details.

We note that Eq. 4 is an approximate expression for the orbital eccentricity.
This together with the above estimation for $L_{\mathrm{circ}}$ gives rise to
a negative result inside the square root in Eq. 4 for a very small population
of stars. As this is a consequence of a breakdown in the above assumption for
that population of stars, we remove them from our sample. Finally, to check
that Eq. refeccentricity gives rise to a reasonable result for the
eccentricity, we took a slightly different approach, (see for instance Eqs. 1,
2 of Mackereth et al., 2019). While the second scheme also fails for stars
with a pericenter less than the softening length of the TNG simulation in
which we do not have a well-defined expression for the gravitational
potential, our results in the overlapping region indicate a satisfactory
correspondence. Owing to this, in what follows, we use Eq. 4 and infer the
eccentricity accordingly.

Figure 11 presents $\beta(r)$ profiles for the entire stellar population (left
column) and three eccentricity-selected sub-samples: $e\leq 0.3$, $0.3\leq
e\leq 0.7$ and $e\geq 0.7$. We note that low-eccentricity stars have much
lower $\beta$ than the parent population. This is because these stars have
nearly-circular orbits, with $\sigma_{r}$ being much lower than the
$\sigma_{\theta}$ and $\sigma_{\phi}$. On the other hand, stars with higher
eccentricities are more on radial orbits than tangential, producing higher
$\beta$ values.

Since the above eccentricity cut led to a substantial variation in the $\beta$
radial profile, to explore its impact further, in Figure 12, we draw the 1D
probability distribution function of eccentricity; from the left to right, we
analyse the distribution for the entirety of stars, for those restricted to
radii $1\leq r/\mathrm{kpc}\leq 10$ and for stars in the range $10\leq
r/\mathrm{kpc}\leq 100$, respectively. From the plot, it is evident that in
most cases, the eccentricity peaks between 0.2-0.4. Furthermore, there is a
signal of bi-modality in the eccentricity distribution. More explicitly, while
in lower distances, the eccentricity is more concentrated on lower values, at
large radii it shows some levels of bi-modality in which we have two peaks in
the eccentricity distribution; one at relatively small radii and another one
in very high values.

Finally, Figure 13 presents the 2D distribution for the radial dependency of
the eccentricity for a sub-sample 4 out of our galaxy sample. In each plot, we
also show the Spearman correlation between the radius and the eccentricity. It
is clearly seen that these two quantities are un-correlated. Furthermore,
different galaxies have their eccentricity peaks at different values and
locations. For example, while galaxies 4 and 25 have a standard eccentricity
peak of 0.2-0.5 in $r\leq 10$ kpc, galaxies 12 and 20 show some levels of bi-
modality in the radial dependency of the eccentricity. It is seen that in
these galaxies there is a very high eccentric distribution at very low radii,
followed by a second peak at $r\leq 10$ kpc. Perhaps such a bi-modal
distribution in the eccentricity is originated by galaxy mergers. However,
this hypothesis must be validated and we leave this to future work.

Figure 13: 2D distribution plot of the radii vs the eccentricity for a sub-
sample four galaxies from our galaxy sample. The Spearman correlation
demonstrates that the eccentricity is uncorrelated with the radii.

###  4.5 The $\beta$ energy profile

Since energy is an integral of motion, the $\beta$-energy profile,
$\beta_{\mathrm{E}}$, may measure the velocity anisotropy in a less biased
manner. It is then intriguing to compare the $\beta$ radial (in this section
referred as $\beta_{\mathrm{r}}$) and energy profiles against each other and
check their similarities.

Figure 14 shows the $\beta_{\mathrm{E}}$ profile for different galaxies in our
sample. To infer $\beta_{\mathrm{E}}$, the stellar energy is read directly
from the TNG outputs. We then generate some energy bins, from the energy min
to its max value, and compute $\beta_{\mathrm{E}}$ for each of these
individually. As the energy range differs from one galaxy to the other, the
energy min and max are not the same for different galaxies. Furthermore, to
ensure that in computing the $\beta_{\mathrm{E}}$, we use almost similar stars
used in computing the $\beta_{\mathrm{r}}$, we imply a spatial cut of $1\leq
r\leq 100$ kpc. Overlaid on each panel, we also present the
$\beta_{\mathrm{r}}$ profile inferred for the median distance in each energy
shell. More explicitly, in every energy shell, used in computing the
$\beta_{\mathrm{E}}$ profile, we read off the median of the stellar distance
to the center and compute the $\beta_{\mathrm{r}}$ at that particular
position. In each panel, the solid-blue line refers to the
$\beta_{\mathrm{E}}$, while the dashed-red line describes the
$\beta_{\mathrm{r}}$ profile. The overall agreement between the
$\beta_{\mathrm{E}}$ and $\beta_{\mathrm{r}}$ is remarkable. There are however
some differences between these profiles that are worth to be quantified. For
this purpose, we write an empirical estimator

$\delta\beta_{\mathcal{ER}}\equiv\frac{\sum_{i}|\beta^{i}_{\mathrm{E}}-\beta^{i}_{\mathrm{r}}|}{\sum_{i}\mid\beta^{i}_{\mathrm{E}}\mid}.$
(5)

Next, in a way similar to what was done in Sec. 4.1.1, we make distinct halo
classes based on the magnitude of $\delta\beta_{\mathcal{ER}}$. We make 3 halo
categories as listed below:

$\bullet$ Category $\mathcal{A_{E}}$: $\delta\beta_{\mathcal{ER}}\leq 0.35$.
The beta radial and energy profiles for halos part of this class are quite
similar. There are 17 halos in this category, including halos
[4,5,6,7,8,9,10,13,14,15,16,17,18,21,22,23,25]. Comparing the halos in this
category with the ones part of $\mathcal{A}$, it is seen in 65% cases, the
halo from class $\mathcal{A}$ is also part of class $\mathcal{A_{E}}$.

$\bullet$ Category $\mathcal{B_{E}}$: $0.35\leq\delta\beta_{\mathcal{ER}}\leq
0.70$. Halos in this class, show a mild difference from the radial to energy
profile. There are 6 halos in this class, including halos [1, 2, 11, 12, 19,
24]. From this list only one case is part of $\mathcal{B}$, one is from
$\mathcal{A}$ while the rest are part of class $\mathcal{C}$.

$\bullet$ Category $\mathcal{C_{E}}$: $\delta\beta_{\mathcal{ER}}\geq 0.70$.
The $\beta_{\mathrm{E}}$ vs $\beta_{\mathrm{r}}$ in this class are not
matched. There are 2 different halos part of this category, including [3, 20]
and both of these halos are part of $\mathcal{C}$.

##  5 Connection to observations

Having presented the $\beta$ radial profile using different cuts, here we make
an in-depth comparison between the above theoretical outcomes and the actual
observational results. There are two different aspects that could be addressed
in our comparison. Firstly, the overall amplitude of the velocity anisotropy
$\beta$ and secondly its radial shape. While the first case is easier to
study, the second aspect is more challenging as it may depend on different
selections/cuts employed for stars which is then less natural. Owing to this,
we make the comparison for different selections, including $\tilde{r}\geq 1$,
$\tilde{r}\geq 1,\tilde{Z}\geq 5$ and $\tilde{r}\geq 1,\tilde{Z}\geq 5$,
[Fe/H] $\leq$ -1.3.

Below, we first point out few different observational studies. We then make a
comparison between our results and different observations.

Figure 14: $\beta$ energy profile vs the $\beta$ radial profile for different
halos in our galaxy sample. For each energy shell, we infer the median
distance of stars and compute the $\beta_{\mathrm{r}}$ for that location. In
inferring the $\beta_{\mathrm{E}}$, we use a spatial cut of $1\leq r\leq 100$
kpc.

###  5.1 Observational constraints on $\beta(r)$

Observations measuring $\beta$ employ various approaches using different
tracers, including K-giants, the blue horizontal-branch (BHB), RR Lyrae, and
subdwarf stars in the solar neighbourhood. Below, we review few different
observational studies.

Sirko et al. (2004) used a sample of blue horizontal-branch (BHB) stars from
the SDSS survey and inferred the radial and tangential velocity dispersions in
the radial distance $5\leq r/\mathrm{kpc}\leq 30$. In Figure 15, we have
converted these values to some estimators for $\beta$ in the aforementioned
interval.

Bond et al. (2010) used a large sample (18.8 million) of main-sequence stars
derived from SDSS and POSS astrometry and estimated the $\beta\simeq 0.65$ in
a radial range $3\leq r/\mathrm{kpc}\leq 13$.

Deason et al. (2012) used distant BHB stars with the Galactocentric distances
in radial range $16<r/\mathrm{kpc}<48$ as the kinematic tracers of the Milky
Way dark halo and infer $\beta\approx 0.5$, in line with local solar
neighborhood studies.

Sohn et al. (2018) used their HST proper motion (PM) measurements from 16
globular clusters (GCs) and find $\beta=0.609^{0.13}_{-0.229}$ in the galactic
distance range from $10-40$ kpc.

Watkins et al. (2019) used a sample of 34 halo GCs from Gaia data spanning in
the radial range 2-21.1 kpc and estimated the velocity anisotropy
$\beta=0.46.^{+0.15}_{-0.19}$.

Lancaster et al. (2019) assembled a very high purity set of 3064 BHBs using
the spectroscopic data from SDSS as well as the astrometric data courtesy from
the Gaia satellite. They computed the $\beta$, see the orange line in their
Figure 9, in the galactic center radii range 7-50 kpc. Their extended data
points allow us to probe the radial evolution of $\beta$. However, since it
only covers 4 data points we present the results in our generic comparison in
Figure 15.

Cunningham et al. (2019) used the data from HALO7D (Halo Assembly in Lambda
Cold Dark Matter: observation in 7 Dimensions) consisting of the KeckII/DEIMOS
spectroscopic and the HST data applied for the main-sequence turnoff stars and
estimate the $\beta$ for the full HALO7D sample as well as its individual
fields, including COSMOS, GOODS-N, GOODS-S and EGS. The details of the data-
points and the $\beta$ values are given in Table 3 of Cunningham et al.
(2019). In Figure 15 we present all of the results for the full-sample as well
as individual ones.

Kafle et al. (2012) used a sample of 4664 BHB stars chosen from SDSS/SEGUE
survey and measured the $\beta$ out to 25 kpc while estimating this up to 60
kpc. There are various trends seen in the data, including a smooth value of
$\beta\sim 0.5$ in the inner part of the halo ($9\leq r/\mathrm{kpc}\leq 12$),
followed by a rather sharp falls off of $\beta$ in the range $r\sim 13-18$ kpc
and a minimum of $\beta=-1.2$ at $r=17$ kpc and an extra rising of $\beta$ at
larger radii. At the outer parts of the halo, in the interval $25\leq
r/\mathrm{kpc}\leq 56$, the actual measurement of $\beta$ is not possible
owing to the lack of the proper motions. They predicted the $\beta\sim 0.5$ in
this interval. In figure 16, top panel, we focus on these data and compare
them with the TNG50 results obtained using different cuts.

King et al. (2015) used a sample of 6174 faint F-type stars from the Hectospec
spectrograph on MMT telescope plus a sample of 3330 BHB stars from SDSS and
estimated the $\beta$ within 6-30 kpc of the MW. The results are presented in
Table 3 of King et al. (2015). In the middle panel of Figure 16, we present
their data and compare them against the simulation results obtained using
different cuts.

Bird et al. (2019) made a sample of 7664 metal-poor K giants using LAMOST
catalog and measured the $\beta$ at the galactocentric radii between 5-100
kpc. As they matched their line-of-sight velocities to the proper motions from
the Gaia data, their final results might be less bias compared with former
literature, such as Kafle et al. (2012); King et al. (2015). In the bottom
panel of Figure 16, we present their results, comparing them against the TNG50
results obtained using different cuts.

Figure 15: Comparison between the $\beta$ radial profile inferred from the
TNG50 simulation, using different selection criteria, and various isolated
data points from different observations. From the left to right, we present
the $\beta$ profile with $\tilde{r}\geq 1$, $\tilde{r}\geq 1,\tilde{Z}\geq 5$
and $\tilde{r}\geq 1,\tilde{Z}\geq 5,$ [Fe/H] $\leq-1.3$, respectively. In
each panel, we also present the data-points from Sirko et al. (2004); Bond et
al. (2010); Deason et al. (2012); Sohn et al. (2018); Watkins et al. (2019);
Lancaster et al. (2019); Cunningham et al. (2019). It is seen that adding the
metallicity cut with [Fe/H] $\leq$ -1.3, somewhat reduces the diversity of
models as it removes the young stars, substantially. On the contrary, the
model diversity is maximal for $\tilde{r}\geq 1$ in which the diversity is
comparable to the error bars from changing the observations. This demonstrates
that including/altering the stellar types may significantly affect the
$\beta$, causing different observations to have different results!

Figure 16: Comparison between the $\beta$ radial profile from the TNG50
simulation with three extended observational data points. Top panel presents
the data from Kafle et al. (2012), the middle panel shows the results from
King et al. (2015) and the bottom one presents the results from Bird et al.
(2019). In each row, we present the TNG50 results using different selection
criteria. It is seen that there are a good overall agreement between the TNG50
results and these extended observations. The model diversity is somewhat
comparable to the error bars in Kafle et al. (2012), while it is smaller than
the ones from King et al. (2015). Finally, while the first two panels show a
prominent dip in the $\beta$ radial profile, the third one does not show any
particular dips. Since proper motion has not been considered in the first two
cases, the third one may be more reliable.

###  5.2 Comparison between the theory and observation

Here we make a comparison between the aforementioned observational results
listed in Sec. 5.1 and the outcome of the TNG50 simulation obtained from
employing different selection criteria.

Before we dig into the details of the comparison, we make some remarks
regarding to the potential impact of substructures on the $\beta$ radial
profile. Observationally, Bird et al. (2019) showed that including the effect
of the Sagittarius stream may lead to a dip in the $\beta$ radial profile.
Consequently, the inclusion and the removal of the substructures are in
principle very important. Theoretically, there have been some studies (see for
e.g Loebman et al., 2018, and references therein) showing that the
substructures may alter the shape of the $\beta$ radial profile. In the
context of the TNG50 simulations, specifically, the substructures are all
removed from the central halos and are only part of the group catalogue. As
the entire focus on this paper is on the central halos, our analysis is safe
of the direct impact of the substructures. Our analysis however is missing the
indirect impact of substructures through their gravitational potential. This
is much harder to quantify as we need to model their orbits which require
greater time cadence in the outputs. We therefore ignore this effect in the
current work and leave its analysis to a future study.

We split the comparison to two different sets; presenting them in two separate
plots.

In Figure 15, we present the individual data points from 7 different
observational data sets (Sirko et al., 2004; Bond et al., 2010; Deason et al.,
2012; Sohn et al., 2018; Watkins et al., 2019; Lancaster et al., 2019;
Cunningham et al., 2019). Adding the TNG50 results obtained from 3 different
selection criteria, $\tilde{r}\geq 1$ (left panel), $\tilde{r}\geq
1,\tilde{Z}\geq 5$ (middle panel) and $\tilde{r}\geq 1,\tilde{Z}\geq 5$,
[Fe/H] $\leq$ -1.3 (right panel). From the plot, it is inferred that the model
diversity is maximal for the left panel, while it is minimal in the right
panel with metallicity cuts. This means that it is easier to find a
theoretical model that matches more observations from the left panel than the
right. Furthermore, despite the fact that there is not any single galaxy that
matches all of different observational points, there are quite few of them
that gets close to quite few observations, especially for the left and middle
panels with no metallicity cuts. To summarize, the diversity of different
realizations from the TNG50 are comparable with the differences between
various observations. This means that TNG50 models can explain different
observations.

In Figure 16, we focus on three extended observations from Kafle et al.
(2012); King et al. (2015); Bird et al. (2019) in which we put the data from
three individual observations in different rows and make the comparison
against the TNG50 results. Again, different panels refer to various selection
criteria. The advantage of using the same observation is that we are generally
less biased by different selection criteria from various observations.
However, we are still affected by some systematic errors owing to the fact
that the observational data of Kafle et al. (2012); King et al. (2015) were
collected before Gaia. This is likely to introduce some biases owing to the
lack of the actual treatment of the proper motion (see e.g. Hattori et al.,
2017, for more details). Consequently, the interpretation of the $\beta$
radial profile might not be straightforward. AS it is seen from the bottom
row, the analysis of Bird et al. (2019) did not lead to a strong dip in the
$\beta$ radial profile.

From the plot, it is seen that there are few models that get very close to
some data points, and within their error bars. Again, the model capabilities
in reproducing the data slightly depend on different selection criteria.

In summary, the TNG50 results are in overall good agreement with different
observational points, though some of these are closer to the data.

Since the diversity of TNG50 models diminish by considering a metallicity cut,
we expect to have more population of models near the peak of agreement while
less of them away from the peak. On the contrary, for models with a higher
diversity, the chance to find a model realization that passes close to the
rare data points are higher, though the population of models near the peak of
agreement is less.

##  6 Conclusions

In this manuscript, we analyzed the stellar kinematics as traced with the
velocity anisotropy, $\beta$, profile for a sample of 25 Milky Way like
galaxies from TNG50 simulation. We made an in-depth study of the $\beta$
radial and the energy profiles, at redshift of $z=0$, exploring the impact of
different selection criteria on the aforementioned quantities.

Below, we summarize the main take-aways from this study:

$\bullet$ It is observed that the stellar metallicity is negatively correlated
with the stellar age, with a Spearman coefficient $|r|\simeq 0.5-0.8$.
Furthermore, it is also negatively correlated with the location of stars,
meaning that older stellar population are mostly located in the galaxy
outskirt while the younger ones are closer to the center. This is expected as
the MW is still in its star-forming stage. Based on this, we explored the
impact of different metallicity cuts on the velocity anisotropy of stars.

$\bullet$ We inferred the $\beta$ radial profile using different selection
criteria and made halo classification, categories $\mathcal{A}$-$\mathcal{C}$,
based on the sensitivity of their profiles on various cuts.

$\bullet$ Halos part of $\mathcal{A}$ class, show a very mild sensitivity to
different selections, where $\delta\beta\leq 0.35$. And members of category
$\mathcal{B}$ exhibit some levels of sensitivity, with
$0.35\leq\delta\beta\leq 0.7$, to different cuts. Finally, halos associated
with $\mathcal{C}$ show a very strong dependencies, $\delta\beta\geq 0.7$, to
model selections. We showed that there are 48%, 16% and 36% of halos in each
of these classes, respectively.

$\bullet$ We analyzed the $\beta$ radial profile and demonstrated that it is
divided to either monotonically increasing profiles (32% of halos) or the ones
with peaks/troughs (68% of population). It is shown that, the $\beta$ profile
from the monotonically increasing class smoothly increases from the interior
to the exterior part of the halo. While, $\beta$ associated with the class
with peaks/troughs experiences some local fluctuations. We further demonstrate
that almost all of the halos from the first class are part of category
$\mathcal{A}$, with only one exception being part of class $\mathcal{B}$. On
the contrary, members of the class with peaks/troughs are mostly part of
category $\mathcal{B}$ or $\mathcal{C}$, with only 29% of them being part of
category $\mathcal{A}$. This means that care should be taken in interpreting
the observed peaks/troughs as they might depend on the actual selection
criteria!

$\bullet$ We probed the impact of different metallicity cuts on the $\beta$
profile and demonstrated that older stars are less rotationally supported and
thus have higher $\beta$. On the other hand, younger stars are more
rotationally supported and have lower $\beta$.

$\bullet$ We explored the impact of different eccentricity based radial cuts
on the $\beta$ radial profile. More explicitly, we divided stars to low,
medium and high eccentric chunks. This explicitly showed that low eccentricity
orbits correspond to lower $\beta$ than the high eccentricity orbits. This is
expected as low eccentric orbits correspond to less radial velocity variance
which naturally leads to smaller $\beta$.

$\bullet$ We computed $\beta$ for different stellar types; demonstrating that
stars on radial orbits have higher velocity anisotropy than those on prograde
and retrograde orbits. This is also expected as radial orbits have lower
absolute angular momentum and thus higher radial velocity variance.

$\bullet$ Moving to the energy space, we calculated the $\beta$ energy profile
and compared that with the $\beta$ radial profile where we made another halo
classification based on their similarities/differences. We explicitly showed
that 65% of halos which present very mild sensitivity on different selection
criteria, in the real space, are among those for which the $\beta$ radial and
energy profiles are very closely related. On the contrary, halos for which the
$\beta$ radial/energy profile are rather different, are entirely parts of
category $\mathcal{C}$. This again establishes that we should be very careful
when drawing any conclusions from the $\beta$ analysis.

$\bullet$ Furthermore, we proposed that one quantitative way to examine
whether the observational results are somewhat robust is to also compute the
$\beta$ energy profile and compare that with the $\beta$ radial profile. Based
on our conclusion, it is likely that those with very similar profiles are more
robust against altering the selection criteria than the ones with different
profiles. However, as already indicated above, the current analysis and
comparison are purely theoretical. From the observational perspective, we need
to estimate the gravitational potential. The current paper is not meant to
show that we identify the unbiased sample of stars using this comparison.
Merely that for a chosen stellar sample, this comparison may indicate how much
the chosen sample might be sensitive to different criteria.

$\bullet$ Ultimately, we overlaid our theoretical outcomes, with few different
selection criteria, on top of some recent observational results and compared
them with each other. Our comparison was made both at the level of isolated
data points, in which each observation gives us only one point, as well as the
extended observations in which the $\beta$ radial profile is reported in more
than one single point for a particular observation. The latter one is more
meaningful since we are not affected by the biases coming from changing the
observations or the selection criteria. It is shown that there are reasonable
agreements between the theoretical predictions and the observational results
both for the isolated data points as well as the extended observations. It is
seen for most of halos there are some data points that easily pass through the
profile.

## Data Availability

Data directly related to this manuscript and its figures is available to share
on reasonable requests from the corresponding author. The IllustrisTNG and
TNG50 simulations are publicly available and accessible at www.tng-
project.org/data (Nelson et al., 2019b).

## acknowledgement

It is a great pleasure to thank Sirio Belli, Sownak Bose, Shep Doeleman and
Sandro Tacchella for the very insightful conversations. We are also grateful
to the referee for very constructive report that improved the quality of this
paper. Razieh Emami acknowledges the support by the Institute for Theory and
Computation at the Center for Astrophysics. We thank the supercomputer
facility at Harvard where most of the simulation work was done. MV
acknowledges support through an MIT RSC award, a Kavli Research Investment
Fund, NASA ATP grant NNX17AG29G, and NSF grants AST-1814053, AST-1814259 and
AST-1909831. FM acknowledges support through the Program "Rita Levi
Montalcini" of the Italian MIUR. The TNG50 simulation was realized with
compute time granted by the Gauss center for Supercomputing (GCS) under GCS
Large-Scale Projects GCS-DWAR on the GCS share of the supercomputer Hazel Hen
at the High Performance Computing Center Stuttgart (HLRS).

Software: matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy
(Oliphant, 2007), seaborn (Waskom et al., 2020), pandas (Reback et al., 2021),
h5py (de Buyl et al., 2016).

## References

  * Abadi et al. (2006) Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2006, MNRAS, 365, 747, doi: 10.1111/j.1365-2966.2005.09789.x
  * Adams et al. (2012) Adams, J. J., Gebhardt, K., Blanc, G. A., et al. 2012, ApJ, 745, 92, doi: 10.1088/0004-637X/745/1/92
  * Adams et al. (2014) Adams, J. J., Simon, J. D., Fabricius, M. H., et al. 2014, ApJ, 789, 63, doi: 10.1088/0004-637X/789/1/63
  * Anders et al. (2017) Anders, F., Chiappini, C., Minchev, I., et al. 2017, A&A, 600, A70, doi: 10.1051/0004-6361/201629363
  * Battaglia et al. (2008) Battaglia, G., Helmi, A., Tolstoy, E., et al. 2008, ApJ, 681, L13, doi: 10.1086/590179
  * Beers et al. (2012) Beers, T. C., Carollo, D., Ivezić, Ž., et al. 2012, ApJ, 746, 34, doi: 10.1088/0004-637X/746/1/34
  * Binney (1980) Binney, J. 1980, MNRAS, 190, 873, doi: 10.1093/mnras/190.4.873
  * Bird et al. (2021) Bird, J. C., Loebman, S. R., Weinberg, D. H., et al. 2021, MNRAS, 503, 1815, doi: 10.1093/mnras/stab289
  * Bird et al. (2019) Bird, S. A., Xue, X.-X., Liu, C., et al. 2019, AJ, 157, 104, doi: 10.3847/1538-3881/aafd2e
  * Bird et al. (2020) —. 2020, arXiv e-prints, arXiv:2005.05980.  https://arxiv.org/abs/2005.05980
  * Bond et al. (2010) Bond, N. A., Ivezić, Ž., Sesar, B., et al. 2010, ApJ, 716, 1, doi: 10.1088/0004-637X/716/1/1
  * Chiba & Yoshii (1998) Chiba, M., & Yoshii, Y. 1998, AJ, 115, 168, doi: 10.1086/300177
  * Cinzano & van der Marel (1994) Cinzano, P., & van der Marel, R. P. 1994, MNRAS, 270, 325, doi: 10.1093/mnras/270.2.325
  * Courteau et al. (2014) Courteau, S., Cappellari, M., de Jong, R. S., et al. 2014, Reviews of Modern Physics, 86, 47, doi: 10.1103/RevModPhys.86.47
  * Cunningham et al. (2016) Cunningham, E. C., Deason, A. J., Guhathakurta, P., et al. 2016, ApJ, 820, 18, doi: 10.3847/0004-637X/820/1/18
  * Cunningham et al. (2019) Cunningham, E. C., Deason, A. J., Sanderson, R. E., et al. 2019, ApJ, 879, 120, doi: 10.3847/1538-4357/ab24cd
  * de Buyl et al. (2016) de Buyl, P., Huang, M.-J., & Deprez, L. 2016, arXiv e-prints, arXiv:1608.04904.  https://arxiv.org/abs/1608.04904
  * Deason et al. (2012) Deason, A. J., Belokurov, V., Evans, N. W., & An, J. 2012, MNRAS, 424, L44, doi: 10.1111/j.1745-3933.2012.01283.x
  * Dejonghe & Merritt (1992) Dejonghe, H., & Merritt, D. 1992, ApJ, 391, 531, doi: 10.1086/171368
  * Dekel et al. (2005) Dekel, A., Stoehr, F., Mamon, G. A., et al. 2005, Nature, 437, 707, doi: 10.1038/nature03970
  * Diakogiannis et al. (2014a) Diakogiannis, F. I., Lewis, G. F., & Ibata, R. A. 2014a, MNRAS, 443, 598, doi: 10.1093/mnras/stu1153
  * Diakogiannis et al. (2014b) —. 2014b, MNRAS, 443, 610, doi: 10.1093/mnras/stu1154
  * Diemand et al. (2005) Diemand, J., Madau, P., & Moore, B. 2005, MNRAS, 364, 367, doi: 10.1111/j.1365-2966.2005.09604.x
  * Dietz et al. (2020) Dietz, S. E., Yoon, J., Beers, T. C., & Placco, V. M. 2020, ApJ, 894, 34, doi: 10.3847/1538-4357/ab7fa4
  * El-Badry et al. (2017) El-Badry, K., Wetzel, A. R., Geha, M., et al. 2017, ApJ, 835, 193, doi: 10.3847/1538-4357/835/2/193
  * Emami et al. (2020) Emami, R., Hernquist, L., Alcock, C., et al. 2020, arXiv e-prints, arXiv:2012.12284.  https://arxiv.org/abs/2012.12284
  * Emami et al. (2021a) Emami, R., Genel, S., Hernquist, L., et al. 2021a, ApJ, 913, 36, doi: 10.3847/1538-4357/abf147
  * Emami et al. (2021b) Emami, R., Hernquist, L., Alcock, C., et al. 2021b, in American Astronomical Society Meeting Abstracts, Vol. 53, American Astronomical Society Meeting Abstracts, 324.07
  * Frankel et al. (2019) Frankel, N., Sanders, J., Rix, H.-W., Ting, Y.-S., & Ness, M. 2019, ApJ, 884, 99, doi: 10.3847/1538-4357/ab4254
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Helmi, A., van Leeuwen, F., et al. 2018, A&A, 616, A12, doi: 10.1051/0004-6361/201832698
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41, doi: 10.1111/j.1365-2966.2005.09321.x
  * Gilmore et al. (2007) Gilmore, G., Wilkinson, M. I., Wyse, R. F. G., et al. 2007, ApJ, 663, 948, doi: 10.1086/518025
  * Grand et al. (2018) Grand, R. J. J., Helly, J., Fattahi, A., et al. 2018, MNRAS, 481, 1726, doi: 10.1093/mnras/sty2403
  * Hani et al. (2019) Hani, M. H., Ellison, S. L., Sparre, M., et al. 2019, MNRAS, 488, 135, doi: 10.1093/mnras/stz1708
  * Hattori et al. (2017) Hattori, K., Valluri, M., Loebman, S. R., & Bell, E. F. 2017, ApJ, 841, 91, doi: 10.3847/1538-4357/aa71aa
  * Hattori et al. (2013) Hattori, K., Yoshii, Y., Beers, T. C., Carollo, D., & Lee, Y. S. 2013, ApJ, 763, L17, doi: 10.1088/2041-8205/763/1/L17
  * Huang et al. (2015) Huang, Y., Liu, X.-W., Zhang, H.-W., et al. 2015, Research in Astronomy and Astrophysics, 15, 1240, doi: 10.1088/1674-4527/15/8/010
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jeans (1915) Jeans, J. H. 1915, MNRAS, 76, 70, doi: 10.1093/mnras/76.2.70
  * Johnson et al. (2021) Johnson, J. W., Weinberg, D. H., Vincenzo, F., et al. 2021, arXiv e-prints, arXiv:2103.09838.  https://arxiv.org/abs/2103.09838
  * Kafle et al. (2012) Kafle, P. R., Sharma, S., Lewis, G. F., & Bland-Hawthorn, J. 2012, ApJ, 761, 98, doi: 10.1088/0004-637X/761/2/98
  * Kewley & Ellison (2008) Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183, doi: 10.1086/587500
  * King et al. (2015) King, Charles, I., Brown, W. R., Geller, M. J., & Kenyon, S. J. 2015, ApJ, 813, 89, doi: 10.1088/0004-637X/813/2/89
  * Kleyna et al. (2001) Kleyna, J. T., Wilkinson, M. I., Evans, N. W., & Gilmore, G. 2001, ApJ, 563, L115, doi: 10.1086/338603
  * Koch et al. (2007) Koch, A., Kleyna, J. T., Wilkinson, M. I., et al. 2007, AJ, 134, 566, doi: 10.1086/519380
  * Krajnović et al. (2005) Krajnović, D., Cappellari, M., Emsellem, E., McDermid, R. M., & de Zeeuw, P. T. 2005, MNRAS, 357, 1113, doi: 10.1111/j.1365-2966.2005.08715.x
  * Lancaster et al. (2019) Lancaster, L., Koposov, S. E., Belokurov, V., Evans, N. W., & Deason, A. J. 2019, MNRAS, 486, 378, doi: 10.1093/mnras/stz853
  * Loebman et al. (2018) Loebman, S. R., Valluri, M., Hattori, K., et al. 2018, ApJ, 853, 196, doi: 10.3847/1538-4357/aaa0d6
  * Łokas (2009) Łokas, E. L. 2009, MNRAS, 394, L102, doi: 10.1111/j.1745-3933.2009.00620.x
  * Łokas et al. (2005) Łokas, E. L., Mamon, G. A., & Prada, F. 2005, MNRAS, 363, 918, doi: 10.1111/j.1365-2966.2005.09497.x
  * Mackereth et al. (2019) Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426, doi: 10.1093/mnras/sty2955
  * Mamon et al. (2013) Mamon, G. A., Biviano, A., & Boué, G. 2013, MNRAS, 429, 3079, doi: 10.1093/mnras/sts565
  * Mamon & Łokas (2005) Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 363, 705, doi: 10.1111/j.1365-2966.2005.09400.x
  * Mashchenko (2015) Mashchenko, S. 2015, arXiv e-prints, arXiv:1504.08273.  https://arxiv.org/abs/1504.08273
  * Merritt (1985) Merritt, D. 1985, AJ, 90, 1027, doi: 10.1086/113810
  * Merritt (2015) —. 2015, ApJ, 814, 57, doi: 10.1088/0004-637X/814/1/57
  * Monachesi et al. (2016) Monachesi, A., Gómez, F. A., Grand, R. J. J., et al. 2016, MNRAS, 459, L46, doi: 10.1093/mnrasl/slw052
  * Monachesi et al. (2019) —. 2019, MNRAS, 485, 2589, doi: 10.1093/mnras/stz538
  * Morrison et al. (1990) Morrison, H. L., Flynn, C., & Freeman, K. C. 1990, AJ, 100, 1191, doi: 10.1086/115587
  * Naidu et al. (2020) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48, doi: 10.3847/1538-4357/abaef4
  * Nelson et al. (2019a) Nelson, D., Pillepich, A., Springel, V., et al. 2019a, MNRAS, 490, 3234, doi: 10.1093/mnras/stz2306
  * Nelson et al. (2019b) Nelson, D., Springel, V., Pillepich, A., et al. 2019b, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science and Engineering, 9, 10, doi: 10.1109/MCSE.2007.58
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196, doi: 10.1093/mnras/stz2338
  * Rashkov et al. (2013) Rashkov, V., Pillepich, A., Deason, A. J., et al. 2013, ApJ, 773, L32, doi: 10.1088/2041-8205/773/2/L32
  * Reback et al. (2021) Reback, J., jbrockmendel, McKinney, W., et al. 2021, pandas-dev/pandas: Pandas 1.3.2, v1.3.2, Zenodo, doi: 10.5281/zenodo.5203279
  * Sales et al. (2007) Sales, L. V., Navarro, J. F., Abadi, M. G., & Steinmetz, M. 2007, MNRAS, 379, 1464, doi: 10.1111/j.1365-2966.2007.12024.x
  * Santucci et al. (2020) Santucci, G., Brough, S., Scott, N., et al. 2020, ApJ, 896, 75, doi: 10.3847/1538-4357/ab92a9
  * Schuster et al. (2012) Schuster, W. J., Moreno, E., Nissen, P. E., & Pichardo, B. 2012, A&A, 538, A21, doi: 10.1051/0004-6361/201118035
  * Sirko et al. (2004) Sirko, E., Goodman, J., Knapp, G. R., et al. 2004, AJ, 127, 914, doi: 10.1086/381486
  * Smith et al. (2009) Smith, M. C., Evans, N. W., Belokurov, V., et al. 2009, MNRAS, 399, 1223, doi: 10.1111/j.1365-2966.2009.15391.x
  * Sohn et al. (2012) Sohn, S. T., Anderson, J., & van der Marel, R. P. 2012, ApJ, 753, 7, doi: 10.1088/0004-637X/753/1/7
  * Sohn et al. (2018) Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862, 52, doi: 10.3847/1538-4357/aacd0b
  * Sommer-Larsen et al. (1994) Sommer-Larsen, J., Flynn, C., & Christensen, P. R. 1994, MNRAS, 271, 94, doi: 10.1093/mnras/271.1.94
  * Stinson et al. (2013) Stinson, G. S., Brook, C., Macciò, A. V., et al. 2013, MNRAS, 428, 129, doi: 10.1093/mnras/sts028
  * Strigari et al. (2007) Strigari, L. E., Bullock, J. S., Kaplinghat, M., et al. 2007, ApJ, 669, 676, doi: 10.1086/521914
  * Thom et al. (2005) Thom, C., Flynn, C., Bessell, M. S., et al. 2005, MNRAS, 360, 354, doi: 10.1111/j.1365-2966.2005.09038.x
  * Tully & Fisher (1977) Tully, R. B., & Fisher, J. R. 1977, A&A, 500, 105
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Walker et al. (2007) Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2007, ApJ, 667, L53, doi: 10.1086/521998
  * Walker et al. (2009) —. 2009, ApJ, 704, 1274, doi: 10.1088/0004-637X/704/2/1274
  * Waskom et al. (2020) Waskom, M., Botvinnik, O., Ostblom, J., et al. 2020, mwaskom/seaborn: v0.10.0 (January 2020), v0.10.0, Zenodo, doi: 10.5281/zenodo.3629446
  * Watkins et al. (2019) Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W. 2019, ApJ, 873, 118, doi: 10.3847/1538-4357/ab089f
  * Wilkinson et al. (2004) Wilkinson, M. I., Kleyna, J. T., Evans, N. W., et al. 2004, ApJ, 611, L21, doi: 10.1086/423619
  * Wolf et al. (2010) Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220, doi: 10.1111/j.1365-2966.2010.16753.x
  * Zuo et al. (2017) Zuo, W., Du, C., Jing, Y., et al. 2017, ApJ, 841, 59, doi: 10.3847/1538-4357/aa70e6