# Early-type galaxy density profiles from IllustrisTNG:
III. Effects on outer kinematic structure

Yunchong Wang1,2,3, Shude Mao1, Mark Vogelsberger4, Volker Springel5, Lars
Hernquist6 and Risa H. Wechsler2,3,7
1Department of Astronomy, Tsinghua University, Beijing, 100084, China
2Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA
94305, USA
3Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450,
Stanford University, Stanford, CA 94305, USA
4Kavli Institute for Astrophysics and Space Research, Department of Physics,
MIT, Cambridge, MA 02139, USA
5Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748,
Garching, Germany
6Institute for Theory and Computation, Harvard-Smithsonian Center for
Astrophysics, 60 Garden Street, Cambridge, MA 02138
7SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA E-mail:
<EMAIL_ADDRESS>
(Accepted 2022 May 13. Received 2022 May 2; in original form 2021 November 17)

###### Abstract

Early-type galaxies (ETGs) possess total density profiles close to isothermal,
which can lead to non-Gaussian line-of-sight velocity dispersion (LOSVD) under
anisotropic stellar orbits. However, recent observations of local ETGs in the
MASSIVE Survey reveal outer kinematic structures at $1.5R_{\mathrm{eff}}$
(effective radius) that are inconsistent with fixed isothermal density
profiles; the authors proposed varying density profiles as an explanation. We
aim to verify this conjecture and understand the influence of stellar assembly
on these kinematic features through mock ETGs in IllustrisTNG. We create mock
Integral-Field-Unit observations to extract projected stellar kinematic
features for 207 ETGs with stellar mass $M_{\ast}\geqslant
10^{11}\mathrm{M_{\astrosun}}$ in TNG100-1. The mock observations reproduce
the key outer ($1.5R_{\mathrm{eff}}$) kinematic structures in the MASSIVE
ETGs, including the puzzling positive correlation between velocity dispersion
profile outer slope $\gamma_{\mathrm{outer}}$ and the kurtosis $h_{4}$’s
gradient. We find that $h_{4}$ is uncorrelated with stellar orbital anisotropy
beyond $R_{\mathrm{eff}}$; instead we find that the variations in
$\gamma_{\mathrm{outer}}$ and outer $h_{4}$ (a good proxy for $h_{4}$
gradient) are both driven by variations of the density profile at the
outskirts across different ETGs. These findings corroborate the proposed
conjecture and rule out velocity anisotropy as the origin of non-Gaussian
outer kinematic structure in ETGs. We also find that the outer kurtosis and
anisotropy correlate with different stellar assembly components, with the
former related to minor mergers or flyby interactions while the latter is
mainly driven by major mergers, suggesting distinct stellar assembly origins
that decorrelates the two quantities.

###### keywords:

galaxies: elliptical and lenticular, cD – galaxies: evolution – galaxies:
structure – galaxies: kinematics and dynamics – methods: numerical

††pubyear: 2021††pagerange: Early-type galaxy density profiles from
IllustrisTNG:  III. Effects on outer kinematic structure–B

##  1 Introduction

Early-type galaxies (ETGs, e.g., Dressler 1980; Djorgovski & Davis 1987) are
recognized as the ‘red and dead’ end products of hierarchical galaxy formation
(Cole et al., 2000; Springel et al., 2001; De Lucia et al., 2006; De Lucia &
Blaizot, 2007). Numerical simulations over the past decade have shown that the
formation path of early-type galaxies is well-represented by a two-phase
scenario (Naab et al., 2007; Guo & White, 2008; Oser et al., 2010; Johansson
et al., 2012; Moster et al., 2013; Rodriguez-Gomez et al., 2016), comprising
an active phase ($z\gtrsim 2$) dominated by gas-rich mergers and bursty in-
situ star formation (Hopkins et al., 2006, 2008a, 2008b; Hopkins et al.,
2009c; Hopkins et al., 2009d; Wellons et al., 2015), as well as a passive
phase ($z\lesssim 2$) dominated by dry mergers and accretion of ex-situ-formed
stars (Nipoti et al., 2009a; Nipoti et al., 2009b; Remus et al., 2013; Wellons
et al., 2016), as well as quenching by Active Galactic Nuclei (AGN) feedback
(Silk & Rees, 1998; King, 2003; Wyithe & Loeb, 2003; Di Matteo et al., 2005;
Springel et al., 2005; Fabian, 2012; Kormendy & Ho, 2013).

An important feature of ETGs found in observations through strong and weak
gravitational lensing (Koopmans et al., 2006; Gavazzi et al., 2007; Koopmans
et al., 2009; Barnabè et al., 2009, 2011; Auger et al., 2010; Ruff et al.,
2011; Sonnenfeld et al., 2013; Li et al., 2018b; Lyskova et al., 2018),
dynamical modeling (Tortora et al., 2014; Cappellari et al., 2015; Serra et
al., 2016; Poci et al., 2017; Bellstedt et al., 2018; Li et al., 2019) and
X-ray observations of gas dynamics (Humphrey et al., 2006; Humphrey & Buote,
2010), is that their total radial matter density profile is well described by
a single power-law model $\rho(r)\propto r^{-\gamma^{\prime}}$ with small
intrinsic scatter around the slope of $\gamma^{\prime}=2$. This feature is
known as the ‘bulge–halo conspiracy’, as neither stars nor dark matter follow
a single power-law model with the slope of $\gamma^{\prime}=2$, but their
combined profile ‘conspires’ to take the form of a Singular Isothermal Sphere
(SIS) (an ideal gas sphere under gravitational-hydrostatic equilibrium):

$\rho(r)=\frac{\sigma^{2}}{2\pi Gr^{2}}\Longleftrightarrow\rho(r)\propto
r^{-\gamma^{\prime}},\gamma^{\prime}=2$ (1)

where $\sigma$ is the one-dimensional velocity dispersion of stars. In the
case of isotropic stellar orbits,
$\sigma=\sigma_{r}=\sigma_{\theta}=\sigma_{\phi}$ where the subscripts denote
the $(r,\theta,\phi)$ directions in spherical coordinates. In this case, a
constant $\sigma$ leads to a constant circular velocity $v_{c}$ with radius:

$v_{c}^{2}(r)=\frac{GM(<r)}{r}=\frac{G}{r}\int_{0}^{r}\rho(r^{\prime})4\pi
r^{\prime 2}dr^{\prime}=2\sigma^{2}.$ (2)

Therefore, the ETG total density profile characterized by $\gamma^{\prime}=2$
is naturally linked to its kinematic structure characterized by a flat
velocity dispersion/circular velocity radial profile under stellar orbit
isotropy. Any deviations from a flat velocity dispersion curve (or the
underlying circular velocity) will indicate a deviation of the total density
profile from an SIS model ($\gamma^{\prime}\neq 2$). If the 3D velocity
dispersion profile can be approximated by a power law $\sigma(r)\propto
r^{\alpha}$, then $\gamma^{\prime}<2$ when $\alpha>0$ and $\gamma^{\prime}>2$
when $\alpha<0$.

However, in the presence of a radial or tangential stellar orbital anisotropy,
the projected stellar velocity dispersion can also vary with radius even when
$\gamma^{\prime}=2$ and $v_{c}$ remains constant with radius. As the
anisotropy cannot be measured directly in observations, neither the 3D
velocity dispersion profile (Binney & Mamon, 1982) nor the density profile
logarithmic slope (Cappellari, 2008; Xu et al., 2017) can be determined in an
unbiased manner without an assumed velocity anisotropy. This causes a
degeneracy in the derived mass profile/velocity dispersion profile and the
assumed velocity anisotropy. Nonetheless, if higher order velocity moments are
measured, which provides non-Gaussian information (especially the kurtosis
$h_{4}$, which is the $4^{\mathrm{th}}$ velocity moment) of the LOSVD, this
degeneracy can in theory be broken through the opposite behavior of $h_{4}$
under radial and tangential velocity anisotropy (Dejonghe & Merritt, 1992;
Gerhard, 1993; Merritt & Saha, 1993; Read & Steger, 2017).

Recent results from the MASSIVE Survey (Ma et al., 2014) found that massive
local ETGs with stellar mass $\gtrsim 10^{11.5}\mathrm{M_{\astrosun}}$ that
have rising velocity dispersion profiles towards their outskirts tend to have
positive $h_{4}$ and positive $h_{4}$ gradient (Veale et al., 2017; Veale et
al., 2018). Ene et al. (2019) also found that most of MASSIVE ETGs have
dropping inner velocity dispersion profiles and increasing $h_{4}$ towards the
galactic center. These trends are in contradiction to the theoretical
expectations that under a fixed total density profile, radial velocity
anisotropy induces a more positive $h_{4}$ accompanied by a decreasing LOSVD
towards the outskirts of a galaxy (vice versa towards the center). To explain
this observed tension, Veale et al. (2018) proposed that the presence of
circular velocity gradients could be the cause, and galaxy-to-galaxy
variations of the total density profile slope ($\gamma^{\prime}$) could be
present in ETGs. Therefore, we aim to verify this conjecture using simulated
ETGs from the state-of-the-art cosmological hydrodynamic simulation
IllustrisTNG (Marinacci et al., 2018; Naiman et al., 2018; Nelson et al.,
2018; Pillepich et al., 2018b; Springel et al., 2018), which possesses a well-
studied ETG sample with near-isothermal total density profiles that are
broadly consistent with their observed counterparts (Paper I, Wang et al.
2020).

If the outer kinematic structures of ETGs are indeed influenced by the
variations in their density profiles, we expect to see a correlation of their
outer kinematic structures with minor mergers. This is due to the fact that
the evolution of the ETG total density profile at $z<0.5$ is mainly driven by
minor mergers as found in Wang et al. (2019) (Paper II, hereafter) and
similarly in earlier works (Johansson et al., 2012; Hilz et al., 2013; Remus
et al., 2013; Sonnenfeld et al., 2014). Interestingly, Greene et al. (2019)
found strong correlations between $h_{4}$ and stellar populations probes (i.e.
metallicity, metallicity gradients) and suggested cumulative minor mergers
might have led to the old-aged, radially-anisotropic ETGs having positive
$h_{4}$ at the outskirts. Since galaxy mergers also affect the velocity
dispersion profile (Bender et al., 1992; Schauer et al., 2014; Nipoti et al.,
2020) and tend to induce radial velocity anisotropy (Romanowsky et al., 2003;
Hilz et al., 2012), we will also investigate the role of minor mergers in the
co-evolution of ETG outer kinematic structure and their density profile using
the merger histories of the simulated ETGs.

This paper is organized as follows: in Section 2 we introduce the simulation
and selection criteria through which we select our ETG catalog, as well as the
methods to mimic observations for extracting kinematic structure information
out of the ETGs; in Section 3 we present the results for the kinematic
properties of our selected ETG sample along with comparisons to observations;
in Section 4 we further explore the physical interpretation for the formation
of outer kinematic structure in ETGs relating to their total density profiles,
stellar assembly histories, and environment; in Section 5 we summarize our
main conclusions and provide an outlook for future directions of work. In the
following analysis, we assume the Planck-2016 flat-$\Lambda$CDM cosmology
(Planck Collaboration et al., 2016) parameters used by the IllustrisTNG
simulations: $h=0.6774$, $\Omega_{\mathrm{m}}=0.3089$,
$\Omega_{\mathrm{b}}=0.0486$, $\Omega_{\mathrm{\Lambda}}=0.6911$, and
$\sigma_{\mathrm{8}}=0.8159$.

##  2 Methodology

###  2.1 IllustrisTNG Simulations

Cosmological simulations have tremendously improved our understanding of
galaxy formation and cosmology over the past few decades (see Vogelsberger et
al. 2020a for a recent review). The Next Generation Illustris
Simulations111www.tng-project.org (Marinacci et al., 2018; Naiman et al.,
2018; Nelson et al., 2018; Pillepich et al., 2018b; Springel et al., 2018),
IllustrisTNG for short, is a recent set of cosmo-magneto hydrodynamic
simulations evolved with the state-of-the-art moving mesh hydrodynamics code
Arepo (Springel, 2010). They advance the merits of the Illustris Simulations
(Vogelsberger et al., 2013; Torrey et al., 2014), and improve the Illustris
models (Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al.,
2015; Nelson et al., 2015) in terms of AGN and stellar feedback physics
(Weinberger et al., 2017; Pillepich et al., 2018a).

The full physics IllustrisTNG simulation suite reproduces many key relations
in observed galaxies, including the galaxy-color bimodality in the Sloan
Digital Sky Survey (Nelson et al., 2018), the fraction of dark matter within
galaxies at $z=0$ (Lovell et al., 2018), the galaxy mass-metallicity relation
(Torrey et al., 2018, 2019) and the intra-cluster metal distribution
(Vogelsberger et al., 2018), the galaxy size-mass relation evolution (Genel et
al., 2018), galaxy morphology transition (Tacchella et al., 2019) and stellar
orbital fraction (Xu et al., 2019), early-type galaxy total density profiles
(Wang et al., 2020), molecular and atomic hydrogen content in low redshift
galaxies (Diemer et al., 2019; Stevens et al., 2019), star formation
activities and quenched fractions (Donnari et al., 2019), ram-pressure
stripping in dense environments (Yun et al., 2019), gas-phase metallicity
gradients in star-forming galaxies (Hemler et al., 2021), as well as AGN
galaxy occupation and X-ray luminosities (Weinberger et al., 2018; Habouzit et
al., 2019; Terrazas et al., 2020). Although some facets of these comparisons
still exhibit discrepancies with observations to different levels, the
significant improvements over Illustris and the multitude of agreement in
galaxy and cluster level properties with observations demonstrates the
predictive power of IllustrisTNG (e.g., predictions of JWST observation for
high redshift galaxies Vogelsberger et al. 2020b). Therefore, we use ETGs
selected from IllustrisTNG to gain insights on the origin of their outer
kinematic structure as seen in the MASSIVE Survey.

The simulation suite of IllustrisTNG comprises 3 cubic boxes with periodic
boundary conditions, i.e. TNG100 (side length $75\,\mathrm{Mpc}/h$, same as
original Illustris), TNG300 (side length $205\,\mathrm{Mpc}/h$), and TNG50
(side length $35\,\mathrm{Mpc}/h$), with overall higher resolution in smaller
boxes and each box contains several runs with different numerical resolutions.
In this paper, we select galaxies from the highest resolution run of the
TNG100 box, which is best for our purpose of studying ETGs with stellar mass
$\gtrsim 10^{11}\mathrm{M_{\astrosun}}$, since it provides a substantial
sample size of ETGs in this mass range with reasonable mass resolution. TNG100
has a baryonic matter mass resolution of $m_{\mathrm{baryon}}=1.4\times
10^{6}\mathrm{M}_{\mathrm{\astrosun}}$ and a dark matter mass resolution of
$m_{\mathrm{DM}}=8.9\times 10^{6}\mathrm{M}_{\mathrm{\astrosun}}$, each with
$1820^{3}$ resolution elements. The softening length scale of dark matter and
stellar particles is $\epsilon=0.74\,\mathrm{kpc}$ (valid for $z\leq 1$,
scales as $(1+z)^{-1}\epsilon$ at $z>1$), whereas the gravitational softening
of the gas cells is adaptive and has a minimum length scale of 0.19 comoving
kpc. The simulation runs of TNG50, TNG100 and TNG300 boxes are now available
for public data access (Nelson et al., 2019).

###  2.2 Sample selection

Galaxies in IllustrisTNG are identified as gravitationally-bound structures
(subhalos) by Subfind (Springel et al., 2001; Dolag et al., 2009) wiin halos
found using the Friends of Friends (FoF) algorithm based on mean particle
separation length. The largest subhalo in a halo together with its baryonic
component is defined as the central galaxy, and all other subhalos in the halo
are defined as satellite galaxies. We select central galaxies in TNG100 with
total stellar masses of $M_{\ast}\geqslant 10^{11}\mathrm{M_{\astrosun}}$
($\gtrsim 7\times 10^{4}$ resolution elements), which covers the stellar mass
range of the MASSIVE ETGs ($M_{\ast}\geq 3\times
10^{11}\mathrm{M_{\astrosun}}$, Veale et al. 2018). Merger trees that trace
the comprehensive assembly histories of galaxies and dark matter halos are
constructed using the algorithm Sublink (Rodriguez-Gomez et al., 2015).

Figure 1: Schematic flow chart of the main steps for our mock-IFU post-
processing pipeline. The input starts from the projected stellar luminosity
map from particle data, and the output is the 2D kinematic structure maps for
each galaxy, i.e. mock IFU data products. The example galaxy (Subfind ID =
257302) shown has typical elliptical geometry in its luminosity map. The
central projected region is well resolved with a large number of high SNR
spaxels, and there are increasing fractions of Voronoi bins towards the
outskirts. The velocity histogram indicates one of the spaxel/Voronoi bin
demonstrating a significant non-Gaussian LOS velocity profile. Comparing the
mean velocity and dispersion maps shows that the galaxy is dominated by
randomized stellar motion, while its skewness $h_{3}$ distribution is anti-
correlated with the mean velocity distribution, and kurtosis $h_{4}$ is
generally positive for most of the spaxels/Voronoi bins in the map.

We follow Paper I and Paper II for galaxy morphology classification, which is
based on the optical luminosity reconstruction approach detailed in Xu et al.
(2017). The optical light of a galaxy is derived using the Stellar Population
Synthesis (SPS) model (Bruzual & Charlot, 2003) based on the metallicity and
age of its stellar particles which are treated as coeval stellar populations
with Chabrier initial mass function (Chabrier, 2003). A projection dependent
dust attenuation is then applied to the galaxy luminosity and we take the SDSS
$r$-band luminosity after dust processing to calculate azimuthally-averaged
galaxy luminosity radial profiles for morphological classification. We fit
single-component models, i.e. the de Vaucouleurs profile
(S$\mathrm{\acute{e}}$rsic $n=4$) and the exponential profile
(S$\mathrm{\acute{e}}$rsic $n=1$), as well as a combined two-component model
of exponential and de Vaucouleurs profiles to the projected radial luminosity
profiles of galaxies. If a galaxy is better-fit (lower minimum $\chi^{2}$) by
the de Vaucouleurs profile, and demonstrates bulge-dominance in the two-
component fit (bulge-to-total ratio $>0.5$) in all three independent
projections (along the $x$, $y$, and $z$ axes of the simulation box), then it
is considered as early-type.

To make our ETG classification more robust according to Integral-Field-Unit
(IFU hereafter) observations (Li et al., 2018a), we fit a single
S$\mathrm{\acute{e}}$rsic profile to the projected luminosity profiles of the
selected ETGs, and keep only those with S$\mathrm{\acute{e}}$rsic indices
satisfying $n_{\mathrm{x}}\geq 2.5$ , $n_{\mathrm{y}}\geq 2.5$, and
$n_{\mathrm{z}}\geq 2.5$ in all three projections simultaneously (Li et al.
2018a selected ETGs based on $n\geq 2.5$). For the MASSIVE ETGs (Ma et al.,
2014), they were selected based on morphology (E and S0 types from Paturel et
al. 2003) without any specific S$\mathrm{\acute{e}}$rsic index cut applied,
and 77 out of 105 galaxies have $n$ ranging from 2 to 6 cross referencing the
NSA catalog. Therefore, we consider our generic ETG selection criteria a
reasonable choice when comparing to MASSIVE ETGs and we arrive at a sample of
221 well-resolved ETGs after the above-mentioned photometric selections. This
will be further reduced to a final sample of 207 ETGs after removing galaxies
that have been through recent major mergers (see Section 3.1).

###  2.3 Mock observations

We summarize the main steps in post-processing simulation particle and catalog
data to retrieve kinematic properties of our mock ETG sample that mimics
kinematic properties from the observational Integral-Field-Unit (IFU)
spectroscopic surveys. Our pipeline is largely based upon the public code
illustris-tools 222https://github.com/HongyuLi2016/illustris-tools (Li et al.,
2016) to make mock IFU observations, with edits for our IllustrisTNG ETG
sample applied. Fig. 1 shows the work flow of our post-processing pipeline.

We find the center of the stellar component of the galaxy by comparing the
center of mass for all stellar particles in the galaxy and its central 20$\%$
of particles (in 3D radial distance). We do this recursively by setting the
central 20$\%$ as the total region considered to be the total region in the
next step, and compare its center of mass to the 20$\%^{2}$ particles center
of mass, until the difference between the two center of masses in a step drops
below 0.01 kpc. This center for the stellar component is within $\sim 0.05$
kpc of the minimum gravitational potential point defined by Subfind for our
ETG sample. We assume the mean velocity of the stellar particles to be the
mass-weighted average of particle velocities within 15 kpc of the stellar
component center found as above. We calibrate coordinates and velocities of
stellar particles with respect to the stellar component center and mean
velocity of their galaxy.

Afterwards, we make 2D projections of the stellar particles in random
directions (along the z axis of the simulation box) to mimic realistic
observational conditions (we have also tried random projections in the x and y
axis, as well as along the edge-on direction of the galaxy, which all largely
preserve the main findings of this paper). We pixelize the SDSS $r$-band
luminosities of the projected stellar particles onto an
$80\,\mathrm{kpc}\times 80\,\mathrm{kpc}$ square aperture centered on the
stellar component center. We set our pixel size to $0.5\,\mathrm{kpc}\times
0.5\,\mathrm{kpc}$, imitating the local ETG population of the MASSIVE survey,
and we place our fiducial ETGs at $z=0.03$ (angular diameter distance 128 Mpc)
for our mock observations. The angular resolution of our mock luminosity map
is $0.806^{\prime\prime}$ at that redshift, while the MASSIVE galaxies were
observed using the Hobby-Eberly Telescope (HET, Hill et al. 2008) at the
McDonald Observatory, where its operational seeing has a FWHM of
$1.5^{\prime\prime}$ ($\sigma_{\mathrm{PSF}}=0.64^{\prime\prime}$) If we
convolve a Gaussian PSF kernel matching the HET seeing, which has a Gaussian
scatter than our mock pixel size, the PSF normalization will boost each
pixel’s value by itself and does not smear out its flux to other neighbouring
pixels. As such, we consider our pixelization has effectively accounted for
the seeing with it being marginally coarser than the expected PSF.

Next, we model the projected luminosity map with the Multi-Gaussian-Expansion
(MGE) formalism (Emsellem et al., 1994) with the publicly available python
package MgeFit333https://www-astro.physics.ox.ac.uk/~mxc/software/
(Cappellari, 2002). The MGE method models the surface brightness of a galaxy
with a stack of (we choose $N=15$) elliptical 2D Gaussians:

$\Sigma(x^{\prime},y^{\prime})=\sum_{k=1}^{N}\frac{L_{k}}{2\pi\sigma_{k}^{2}q_{k}^{\prime}}\mathrm{exp}\left[-\frac{1}{2\sigma_{k}^{2}}\left(x^{\prime
2}+\frac{y^{\prime 2}}{q_{k}^{\prime 2}}\right)\right],$ (3)

where $L_{k}$, $\sigma_{k}$, $q_{k}^{\prime}$ are the normalization, standard
deviation, and axis ratio of the $k$-th Gaussian. The primed coordinates
$x^{\prime}$, $y^{\prime}$ are the projected 2D angular coordinates (placing
the galaxy at 128 Mpc) in a system where the origin is at the center of the
galaxy and $x^{\prime}$ axis being aligned with the galaxy’s major axis
(Cappellari, 2002). The center and major (minor) axis of the galaxy’s
projected luminosity map are calculated using the top $10\%$-brightest pixels
(Li et al., 2016), with the center being determined by luminosity-weighting
and the major (minor) axis following the eigenvectors of the 2D inertial
tensor (Allgood et al., 2006) spanned by the selected pixels.

The luminosity profile is then sampled in equal angular bins and logarithmic
radial bins to perform the MGE fit. The result of this step produces an
analytical best-fit description of the galaxy surface brightness map with 2D
Gaussians, and also provides realistic effective radius (half-light radius) of
our mock ETGs that mimics observations by integrating the best-fit surface
brightness profiles. A comparison between the MGE-derived effective radius
$R_{\mathrm{eff}}^{\mathrm{MGE}}$ and the effective radius obtained by
directly projecting all stellar particles assigned to the galaxy by Subfind
($R_{\mathrm{eff}}$) is shown in Appendix A, which shows that
$R_{\mathrm{eff}}^{\mathrm{MGE}}$ gives a more realistic description of galaxy
sizes removing much of the intra-cluster light especially in more massive
ETGs. This choice of effective radius is also more flexible to capture
different distinct components in the ETG luminosity profile compared to
simpler half-light-ellipse or single S$\mathrm{\acute{e}}$rsic fits often
adopted in observations (see Section 3.2 in Ma et al. 2014). In the rest of
this paper, we refer to $R_{\mathrm{eff}}^{\mathrm{MGE}}$ when we mention the
effective radius of our simulated ETGs.

Figure 2: The distribution of stellar mass (left panel) and effective radius
(right panel) versus the velocity dispersion measured within the effective
radius. The open markers indicate those ETGs that are either undergoing major
mergers found via their kinematic maps, or that have been through a recent one
(mass ratio $>1/3$) in the previous snapshot recorded in the merger tree. We
remove them in the following analysis and figures following Veale et al.
(2018), and only include the filled marker points in our samples below (207
galaxies remaining out of 221). The histograms at the sides in each subplot
show the marginalized distributions of these three quantities (histograms
include the major merger galaxies). Dashed lines in both panels indicate
scaling relations (fixed slope and fitted intercept) for ETGs as expected from
observations. ‘Plus’ markers in both panels indicate the MASSIVE sample in
Veale et al. (2018). Our simulated ETG sample has only 18 galaxies with
$\log_{10}(M_{\ast}/\mathrm{M_{\astrosun}})\geqslant 11.6$) versus 89 in
MASSIVE, owing to the volume limit of the simulation. The effective radius for
MASSIVE ETGs quoted here are from the 2MASS catalog (which underestimates
galaxy sizes by a factor $\sim 2\ (0.3\ \mathrm{dex})$ compared to the NSA
catalog, but covers the entire MASSIVE sample) as summarized in Ma et al.
(2014). We also show a dotted line in the right panel that indicates a scaling
relation with half the effect radius of the fitted IllustrisTNG ETG scaling
which coincides well with the MASSIVE galaxies.

Furthermore, we generate mock IFU maps for our simulated ETGs focused on their
central regions satisfying:

$r^{\prime}=\sqrt{x^{\prime 2}+\frac{y^{\prime
2}}{q_{\mathrm{p}}^{2}}}<\left(\frac{R_{\mathrm{max}}}{\sqrt{q_{\mathrm{p}}}}\right),$
(4)

where
$R_{\mathrm{max}}=\mathrm{max}\\{2.5R_{\mathrm{eff}}^{\mathrm{MGE}},20\,\mathrm{kpc}\\}$
is our mock IFU aperture that mimics observations (Veale et al., 2017). This
aperture setting guarantees that the kinematic maps we generate sample the
outskirts of the galaxies well. After we select these central pixels, we
assume that each individual pixel corresponds to an IFU spaxel (single fiber
in an IFU bundle). We use the Convex-Hull method to efficiently identify
stellar particles that are projected within our IFU aperture, and assign them
to their nearest spaxel by querying the KD-Tree constructed for all spaxel
anchor points. Following that, we Voronoi-bin the outskirt spaxels to the
target signal-to-noise ratio (SNR) of 20 to ensure kinematic information
quality, and leave single spaxels in the central region with high SNR
unbinned.

Finally, we calculate the line-of-sight (LOS) stellar velocity distribution
for each spaxel/Voronoi-bin in the 2D IFU aperture. If a spaxel/Voronoi-bin
has $N$ stellar particles projected within it, then we construct the LOS
velocity histogram of these particles using $\left[\sqrt{N}\right]$ equi-
velocity bins. We perform a least-squares fit to the stellar LOS velocity
histogram in each spaxel/Voronoi-bin with a fourth-order Gauss-Hermite
function to extract the mean ($\bar{v}$), dispersion ($\sigma$), skewness
($h_{3}$), and kurtosis ($h_{4}$) of the LOS velocity distribution:

$\begin{split}f(v)&=\frac{1}{\sqrt{2\pi\sigma^{2}}}\,\mathrm{exp}\left(-\frac{x^{2}}{2}\right)\\\
&\quad\quad\times\left[1+h_{3}\mathrm{H_{3}}\left(x\right)+h_{4}\mathrm{H_{4}}\left(x\right)\right],x=\frac{v-\bar{v}}{\sigma},\end{split}$
(5)

where $\bar{\nu}$ and $\sigma$ are the mean velocity and dispersion in the
spaxel/Voronoi bin. $\mathrm{H_{3}}(x)=(2/\sqrt{3})x^{3}-\sqrt{3}x$ and
$\mathrm{H_{4}}(x)=(\sqrt{6}/3)x^{4}-\sqrt{6}x^{2}+\sqrt{6}/4$ are the third
and fourth order normalized Gauss-Hermite functions, while their coefficients
contain the important kinematic information of the third velocity moment
$h_{3}$ (skewness) and the fourth velocity moment $h_{4}$ (kurotsis).
Qualitatively, a positive $h_{3}$ indicates asymmetric bias of the velocity
distribution towards velocities less than the mean, and vice versa for
negative $h_{3}$. A positive $h_{4}$ indicates symmetric deviation from a
normal Gaussian with larger tails and a more centrally peaked velocity
distribution, while negative $h_{4}$ indicates smaller tails and centrally
flat (sometimes double peaked) velocity distribution (van der Marel & Franx,
1993). To wrap up, we store the 2D map of the mean, dispersion, skewness, and
kurtosis of the line-of-sight (LOS) velocity distribution for each individual
ETG, and proceed to obtaining radial profiles of these kinematic properties in
the next section.

##  3 Results of kinematic features

Figure 3: From top to bottom, the profiles of the line of sight mean velocity
($\langle v_{\mathrm{LOS}}\rangle$), velocity dispersion
$\sigma_{\mathrm{LOS}}$, skewness ($h_{3}$), and kurtosis ($h_{4}$) are shown
for each individual galaxy in our ETG sample. Each curve represents one
galaxy’s flux-weighted radial profile of the four kinematic quantities. The
left row shows the profile with respect to the projected physical distance to
the galaxy center, while the right row shows the results w.r.t. the projected
distance in units of the effective radius $R_{\mathrm{eff}}^{\mathrm{MGE}}$
for each individual galaxy. The vertical dotted lines in the left column
indicate the location of $R=2$ kpc and $R=20$ kpc, while the vertical dotted
lines in the right column indicate the location of
$R=2.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. We select three individual galaxies
that have rising (solid), flat (dotted-dashed), dropping (dashed) outer
velocity dispersion profiles at 20 kpc (see Section 3.2 for definitions), and
highlight their kinematic profiles in each panel.

###  3.1 The galaxy sample

In this section, we present the fundamental properties of our selected ETG
sample. Fig 2 shows the stellar mass ($M_{\ast}$), the effective radius
derived with the best-fit MGE model ($R_{\mathrm{eff}}^{\mathrm{MGE}}$), and
the mass-weighted average LOS velocity dispersion
($\sigma_{\mathrm{eff}}^{\mathrm{MGE}}$) within
$R_{\mathrm{eff}}^{\mathrm{MGE}}$. The sample is a typical massive (median
stellar mass $\log_{10}\left(M_{\ast}/\mathrm{M_{\astrosun}}\right)=11.18$)
ETG sample with well-resolved spatial extent (median effective radius of 7
kpc), and dynamically hot with significant velocity dispersion (median
velocity dispersion $162\,\mathrm{km\,s^{-1}}$). The LOS velocity dispersion
within the effective radius positively correlates with stellar mass following
the Tully–Fisher relation $M_{\ast}\propto\sigma^{4}$ (McGaugh et al., 2000;
Bell & de Jong, 2001; Williams et al., 2010). It also correlates with the
effective radius, which is observed to roughly follow $R_{\mathrm{eff}}\propto
M_{\ast}^{0.5}$ in the local universe (Cassata et al., 2013; Huertas-Company
et al., 2013; Xie et al., 2015), and leads to a scaling with the LOS stellar
velocity dispersion in the form of $R_{\mathrm{eff}}\propto\sigma^{2}$. We fit
these two scaling relations to our IllustrisTNG ETG sample with fixed power-
law slopes and free intercepts (shown as dashed lines in Fig. 2). We also plot
in Fig. 2 the stellar masses, effective radii, and the average LOS stellar
velocity dispersion for MASSIVE galaxies in Veale et al. (2018). The
$M_{\ast}$-$\sigma$ scaling of MASSIVE galaxies is rather consistent with the
IllustrisTNG ETGs, although IllustrisTNG approaches the volume limit for
massive galaxies (only 18 galaxies with
$\log_{10}(M_{\ast}/\mathrm{M_{\astrosun}})\geqslant 11.6$) versus the MASSIVE
Survey (89 galaxies with $\log_{10}(M_{\ast}/\mathrm{M_{\astrosun}})\geqslant
11.6$). The effective sizes quoted in the right panel of Fig. 2 are from the
2MASS catalog which covers essentially all MASSIVE ETGs. However, due to the
shallower survey depth of 2MASS, galaxy sizes are underestimated by a factor
$\sim 2$ due to insufficient sensitivity to galaxy outskirts compared to sizes
from the NSA catalog which is based on SDSS DR8 photometry (see Section 3.2 in
Ma et al. 2014). We also plot a scaling relation that shifts the IllustrisTNG
ETG $R_{\mathrm{eff}}$-$\sigma$ scaling by half the effective radii in Fig. 2
which coincides well with the MASSIVE ETGs.

Before diving in to the results, we would like to point out that the selection
of galaxy samples in the MASSIVE Survey (Ma et al., 2014; Veale et al., 2018)
has deliberately removed galaxies in ongoing mergers or show complex merger
remnant structures. We therefore also remove such out-of-equilibrium galaxies
from our sample to further make fair comparisons with the observations. First,
we remove galaxies at $z=0$ that are currently going through major mergers
apparent from their kinematic structure (mainly LOS mean velocity) and
luminosity maps. Second, we also trace every ETG in our sample back one
snapshot from $z=0$ to $z\approx 0.01$ (140 Myrs ago) along their main
progenitor branch, and remove all galaxies that are remnants of recent major
mergers in that snapshot with merger stellar mass ratio $\mu_{\ast}\geqslant
1/3$. Combining these two removal criteria, our original 221 ETG sample
reduces to 207 galaxies, and the removed galaxies are indicated by empty
diamonds in Fig. 2. We proceed with our analysis in the following with the
final 207 galaxies without significant perturbation.

Figure 4: The distribution of the inner and outer slopes of the velocity
distribution profile. The point (0, 0) indicates the case for a single-power-
law isothermal density profile. The colored points indicate simulated ETGs
with rising (blue), flat (green), and dropping (red) outer velocity dispersion
profile slopes ($\gamma_{\mathrm{outer}}$). The error bars of the points are
the combined fitting errors of the velocity dispersion profile for the inner
and outer slopes. Observation values from the MASSIVE Survey (Veale et al.,
2018) are shown for comparison with plus sign markers.

In Fig. 3 we show the projected profiles of the mean velocity ($\langle
v_{\mathrm{LOS}}\rangle$), the stellar velocity dispersion
($\sigma_{\mathrm{LOS}}$), the skewness ($h_{3}$), and kurtosis ($h_{4}$) of
our 207 ETGs sample. We bin the spaxels/Voronoi bins of the projected map
radially for each galaxy in 20 equal-number bins, and calculate the flux-
weighted (SDSS $r$-band) value of the four quantities in each bin to represent
the value for the average velocity profile in that radial bin. The left column
has the projected physical distance $R$ as the x axis, while the right column
has the scaled projected distance in units of
$R_{\mathrm{eff}}^{\mathrm{MGE}}$ as x-axis. We can see that in both radial
scalings the mean velocity and the skewness profiles are largely flat, while
the LOS velocity dispersion profile takes the form of a variety of shapes, and
that the $h_{4}$ profiles show a mildly increasing trend towards larger radii.

###  3.2 Inner and outer slopes of the velocity dispersion profile

To measure the slopes of the projected velocity dispersion profile and compare
with the results of Veale et al. (2018), we follow their modeling by fitting a
double-power-law model to the dispersion profile. The dispersion profile model
takes the form:

$\sigma_{\mathrm{fit}}(R)=\sigma_{0}2^{\gamma_{1}-\gamma_{2}}\left(\frac{R}{R_{\mathrm{B}}}\right)^{\gamma_{1}}\left(1+\frac{R}{R_{\mathrm{B}}}\right)^{\gamma_{2}-\gamma_{1}},$
(6)

where $\gamma_{1}$ and $\gamma_{2}$ are the asymptotic slopes of the profile
towards zero and infinity, $R_{\mathrm{B}}$ the break radius is set to 5 kpc
following Veale et al. (2018), and $\sigma_{0}$ is the normalization of the
dispersion profile. We have verified that the break radius of 5 kpc is a
reasonable description of the IllustrisTNG ETGs. We fit each galaxy’s radial
velocity dispersion profile in the radial range [$2\,\mathrm{kpc}$,
$R_{\mathrm{max}}$] (see Equation 4 for $R_{\mathrm{max}}$) to avoid core
softening in the central regions of the ETGs due to the resolution limit of
the simulation. With an analytic double-power-law description of the
dispersion profiles, we can define the logarithmic slope of the fitted
dispersion profile at any radius $R$ as:

$\gamma(R)=\frac{d\mathrm{log}\,\sigma_{\mathrm{fit}}(R)}{d\mathrm{log}R}=\gamma_{1}+\left(\gamma_{2}-\gamma_{1}\right)\left(\frac{R}{R_{\mathrm{B}}}\right)\left(1+\frac{R}{R_{\mathrm{B}}}\right)^{-1}.$
(7)

Since in the MASSIVE Survey, the inner and outer slopes of the dispersion
profile are measured at $R=2$ kpc and $R=20$ kpc, we follow that definition
and obtain the inner ($\gamma_{\mathrm{inner}}$) and outer
($\gamma_{\mathrm{outer}}$) slopes of the velocity dispersion profile:

$\gamma_{\mathrm{inner}}=\gamma(2)=\frac{5}{7}\gamma_{1}+\frac{2}{7}\gamma_{2},\
\gamma_{\mathrm{outer}}=\gamma(20)=\frac{1}{5}\gamma_{1}+\frac{4}{5}\gamma_{2}.$
(8)

Using the slopes $\gamma_{1}$ and $\gamma_{2}$ obtained from the best double-
power-law fit in Equation 6, we can derive the $\gamma_{\mathrm{inner}}$ and
$\gamma_{\mathrm{outer}}$ values for our ETGs accordingly. We also follow
uncertainties introduced in the fitting procedure to $\gamma_{1}$ and
$\gamma_{2}$, and propagate them appropriately to uncertainties in
$\gamma_{\mathrm{inner}}$ and $\gamma_{\mathrm{outer}}$. Galaxies consistent
with $\gamma_{\mathrm{outer}}=0$ within the uncertainties are defined having
‘flat’ outer dispersion profiles (0 within upper and lower bounds of
$\gamma_{\mathrm{outer}}$), galaxies with lower bound of
$\gamma_{\mathrm{outer}}>0$ are defined as having ‘rising’ outer dispersion
profiles, and galaxies with upper bound $\gamma_{\mathrm{outer}}<0$ are
defined as having ‘dropping’ outer dispersion profiles, consistent with Veale
et al. (2018).

The distribution of the inner and outer slopes (along with their
uncertainties) of the dispersion profiles of our ETG sample are shown in Fig.
4. Galaxies with flat, rising, or dropping outer dispersion profiles are
indicated with differently-colored markers, and we also over-plot the best-fit
values of the observed MASSIVE ETGs from Veale et al. (2018). From the figure
we see that the $\gamma_{\mathrm{outer}}$ of the IllustrisTNG ETGs generally
agree with the range of outer slopes in the observations. The simulation also
has similar levels of uncertainties in the slopes ($\sim 0.1$, not shown for
the observed values) as the observations. The range of
$\gamma_{\mathrm{inner}}<0$ values for our ETGs is also very similar to
observations, however, it is noticeable that there are apparently more
galaxies having $\gamma_{\mathrm{inner}}>0$ within the uncertainties in our
sample, which could be produced by systematic underestimation of the central
stellar velocity dispersion. Since we mainly focus on the behavior of
$\gamma_{\mathrm{outer}}$ in this work, the discrepancy in
$\gamma_{\mathrm{inner}}$ with observations does not affect our findings
below. Nevertheless, this systematic underestimation of the central stellar
velocity dispersion again points to the fact that IllustrisTNG galaxies are
subject to higher central dark matter fractions (Lovell et al., 2018) and halo
contraction (Wang et al., 2020), which adds another galaxy property that needs
to be anchored by future improvements of the underlying feedback physics model
of the simulation.

###  3.3 Higher order velocity moments $h_{3}$ and $h_{4}$

Figure 5: The average kurtosis (left panel, measured within the effective
radius) and kurtosis gradient (right panel, measured within
$0.5R_{\mathrm{eff}}^{\mathrm{MGE}}<R<1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$)
versus the velocity dispersion profile outer slope $\gamma_{\mathrm{outer}}$.
The colored dots indicate the three types of simulated ETGs with different
outer slopes of their velocity dispersion profiles, while the open markers are
observations from the MASSIVE Survey. The point (0, 0) in the left panel
indicates the point where the line-of-sight velocity distribution is a perfect
Gaussian and independent of the projected radius. In the right panel, the
dotted line indicates the best linear fit for the MASSIVE ETGs (from Fig. 6 in
Veale et al. 2018, slope $b=0.096\pm 0.037$, $p=0.012$), and the blue dotted-
dashed line indicates the best linear fit for the IllustrisTNG ETGs (slope
$b=0.092\pm 0.035$, $p=0.009$, shaded region stands for $68\%$ confidence
interval). Clearly, the IllustrisTNG ETGs also show positive correlation
between the kurtosis gradient and outer slope of the velocity dispersion
profile, but there is larger scatter in the $h_{4}$ gradient compared to
observations. The TNG ETGs extend not as far as their MASSIVE Survey
counterparts in the $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ relation due to the
less abundant sampling of the galaxy (Fig. 3) and halo (Fig. 9) masses in
TNG100. Figure 6: Upper panel: The correlation between the skewness $h_{3}$
and the rotation to dispersion ratio $v/\sigma$ at the pixel level measured
around three different radii. Lower panel: The correlation between the average
kurtosis $h_{4}$ and the spherical velocity anisotropy parameter
$\beta_{\mathrm{3D}}$ measured in 3D apertures of different radii. Solid
curves indicate the mean while the shaded regions/error bars represent the
$1\sigma$ confidence intervals in both panels.

In this section, we make direct comparison of the results of $h_{4}$ and
$\gamma_{\mathrm{outer}}$ for IllustrisTNG ETGs with Figure 6 in Veale et al.
(2018), which is the pivotal comparison for this work. With the $h_{4}$
profiles for individual galaxies derived in the previous section, we compute
the average kurtosis within the effective radii $\langle
h_{4}\rangle_{R_{\mathrm{eff}}^{\mathrm{MGE}}}$ with the SDSS $r$-band
luminosity weighted value of $h_{4}$ in all spaxels/Voronoi bins projected
within $R_{\mathrm{eff}}^{\mathrm{MGE}}$ of each galaxy. To derive the
gradient of the kurtosis $\Delta h_{4}(R)/\Delta\mathrm{log}R$, we adopt the
definitions in the observations (Veale et al., 2017; Veale et al., 2018) and
define:

$\frac{\Delta
h_{4}(R)}{\Delta\mathrm{log}R}=\frac{h_{4}(1.5R_{\mathrm{eff}}^{\mathrm{MGE}})-h_{4}(0.5R_{\mathrm{eff}}^{\mathrm{MGE}})}{\mathrm{log_{10}}(1.5R_{\mathrm{eff}}^{\mathrm{MGE}})-\mathrm{log_{10}}(0.5R_{\mathrm{eff}}^{\mathrm{MGE}})},$
(9)

which is the luminosity-weighted average slope of $h_{4}$ within the
logarithmic radial range of [$0.5R_{\mathrm{eff}}^{\mathrm{MGE}}$,
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$]. The spaxels/Voronoi bins included in
the calculation of $h_{4}(R)$ at $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ are
taken from a thin annulus $\pm 0.5$ kpc that of $R$. Due to limited amount of
spaxels available near the galactic center, we compute $h_{4}(R)$ at
$0.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ as the luminosity-weighted average
$h_{4}$ for all spaxels within $R\leqslant
0.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. The results of the mean kurtosis $\langle
h_{4}\rangle_{R_{\mathrm{eff}}^{\mathrm{MGE}}}$ and the kurtosis gradient
$\Delta h_{4}(R)/\Delta\mathrm{log}R$ versus the velocity dispersion profile
outer slope $\gamma_{\mathrm{outer}}$ for our ETG sample are shown in Fig. 5.
In both panels, observation values from Veale et al. (2018) are over-plotted
for comparison to the simulation results. In the right panel, we also show the
best-linear-fit comparison between the IllustrisTNG and MASSIVE $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ relations.

In the left panel of Fig. 5, the IllustrisTNG ETGs’ average kurtosis are
consistent with the range of $h_{4}$ values from observations. In the right
panel of Fig. 5, although the kurtosis gradients of the simulated ETGs have
larger scatter than the observed ETGs, it is remarkable that the simulation
reproduces the key positive correlation between $\Delta
h_{4}(R)/\Delta\mathrm{log}R$ and $\gamma_{\mathrm{outer}}$ as in Figure 6 of
Veale et al. (2018). The dotted-dashed blue line represents the best linear
fit to the IllustrisTNG ETGs (slope $b=0.092\pm 0.035$, $p=0.009$, shaded
region stands for $1\sigma$ confidence interval), and the dotted black line
indicate the best linear fit to the MASSIVE ETGs (slope $b=0.096\pm 0.037$,
$p=0.012$, see Fig. 6 in Veale et al. 2018). Thus, both the slope and
correlation significance for IllustrisTNG ETGs are statistically consistent
with the MASSIVE ETGs, which justifies the use of our mock ETG sample as a
benchmark for understanding the origin of $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ correlation in the
outer kinematic structure of ETGs.

To further check the consistency of higher order kinematic properties of our
simulated ETGs with observations (e.g., Figure 10 in Veale et al. 2017), we
show the relation of the skewness $h_{3}$ versus the rotation to dispersion
ratio $\langle v_{\mathrm{LOS}}\rangle/\sigma_{\mathrm{LOS}}$ in the upper
panel of Fig. 6. Values of $h_{3}$ and $v/\sigma$ for individual
spaxels/Voronoi bins from the different radial apertures (scaled to
$R_{\mathrm{eff}}^{\mathrm{MGE}}$) in different galaxies are stacked together.
Unlike $h_{4}$ which is the fourth velocity moment and has even parity, the
reason why we keep track of the individual spaxel/Voronoi bin level $h_{3}$
information is due to its odd parity, otherwise projection in circular
apertures would average out any underlying signal. The anti-correlation
between $h_{3}$ and $v/\sigma$ seen for the simulated ETGs is consistent with
the theoretical expectation that galaxies with substantial rotation (large
$|v/\sigma|$) will have a large low velocity tail due to the projection effect
of distant background and foreground stars along the line-of-sight, especially
when the galaxy is viewed edge-on. The magnitude of $|h_{3}|\sim 0.1$ is
prevalent from small ($\sim 1R_{\mathrm{eff}}^{\mathrm{MGE}}$) to large ($\sim
2R_{\mathrm{eff}}^{\mathrm{MGE}}$) radii, which also shows consistency with
observations (Krajnović et al., 2011; Veale et al., 2017). Thus, both $h_{3}$
and $h_{4}$ measurements of the IllustrisTNG ETG sample from our mock IFU
pipeline resembles that of the observed MASSIVE sample, indicating that
physical interpretations (Section 4) of these kinematic structures analyzed
from the simulation can be confidently applied to understanding the formation
process of these features in real world ETGs.

###  3.4 De-correlated outer kurtosis and velocity anistropy

The core puzzle in Veale et al. (2018) that sparked the current study is the
positive correlation between $\Delta h_{4}(R)/\Delta\mathrm{log}R$ and
$\gamma_{\mathrm{outer}}$. If an intrinsically flat velocity dispersion
profile (assuming no ordered rotation in an ideal ETG, $v/\sigma$=0) is
affected by the projection effects of velocity anisotropy at the outskirts, we
expect to find radial anisotropy driving positive $h_{4}$ and ‘dropping’
$\gamma_{\mathrm{outer}}$ (vice versa for the tangential case), producing an
anti-correlation of $\Delta h_{4}(R)/\Delta\mathrm{log}R$ and
$\gamma_{\mathrm{outer}}$ (Gerhard et al., 1998). This theoretical expectation
is exactly opposite to what was observed in Veale et al. (2018). A crucial
assumption of this prediction is that they were evaluated under fixed density
profiles (Gerhard, 1993), typically comprising a singular isothermal sphere
halo potential (fixed $\gamma^{\prime}=2$) and a spherical stellar
distribution that have a steeper density profile (fixed $\gamma^{\prime}=3$ or
$4$), which naturally leads to Veale et al. (2018) proposing variations in the
total density profile as an explanation.

Before presenting the connection between the outer kinematic structure and the
density profile, we show in the lower panel of Fig. 6 how the velocity
anisotropy ($\beta_{\mathrm{3D}}$) correlate with $h_{4}$ in our ETGs. The
spherical velocity anisotropy is defined as (Binney & Tremaine, 2008):

$\beta_{\mathrm{3D}}=1-\frac{\sigma_{\mathrm{\phi}}^{2}+\sigma_{\mathrm{\theta}}^{2}}{2\sigma_{\mathrm{r}}^{2}}\,,$
(10)

where $\sigma_{\mathrm{\phi}}$, $\sigma_{\mathrm{\theta}}$, and
$\sigma_{\mathrm{r}}$ are the velocity dispersion in the azimuthal, polar, and
radial directions. As shown in Fig. 6, there is almost no correlation between
$\beta_{\mathrm{3D}}$ and $h_{4}$ averaged in across different aperture sizes.
Particularly, we find a positive $h_{4}$ (median $\sim 0.04$) even for
tangential anisotropy ($\beta<0$), and only a slight positive trend consistent
across different radii at $\langle\beta_{\mathrm{3D}}\rangle_{R}\gtrsim 0.1$.
This suggests that positive $h_{4}$ towards the galactic outskirts is unlikely
driven by systematic radial velocity anisotropy prevalent in ETGs.

In concordance with the conjecture in Veale et al. (2018), we elaborate in
Section 4.1 that the positive $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ correlation is in fact
driven by the variation of their total density profiles across our ETG sample.
The variations in the outer kurtosis and velocity anisotropy also show hints
to be driven by different routes of stellar assembly (Section 4.2). The de-
correlated behavior between the outer $h_{4}$ and $\beta_{\mathrm{3D}}$
($R\geqslant R_{\mathrm{eff}}^{\mathrm{MGE}}$) shown in Fig. 6 can then be
accounted for by the dominant effect from gradients in circular velocity
(Gerhard, 1993; Baes et al., 2005), where the effect of varying density
profiles can overwhelm and randomize the rather weak correlation between
$h_{4}$ and velocity anisotropy.

##  4 Linking kinematic structure with density profile, assembly history and
environment

As shown in the previous section, the IllustrisTNG ETGs have consistent outer
slopes of the velocity dispersion profile ($\gamma_{\mathrm{outer}}$) and
higher order velocity moments ($h_{3}$, $h_{4}$) compared to observations. The
simulated ETGs also reproduce the observed positive correlation between
$\Delta h_{4}(R)/\Delta\mathrm{log}R$ and $\gamma_{\mathrm{outer}}$, and that
velocity anisotropy is insufficient to explain the formation of this trend. In
this section, we explore the validity of circular velocity gradients as a
driving factor of this trend of following the hypothesis in Veale et al.
(2018). We also explore the connection of their kinematic structure to their
stellar assembly and environment.

Figure 7: Upper left: the relation between the $h_{4}$ gradient and the outer
$h_{4}$ values measure at around $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. Upper
right: the relation between the outer $h_{4}$ values and the outer velocity
dispersion profile slopes $\gamma_{\mathrm{outer}}$. Lower left: the relation
between outer $h_{4}$ and the mass-weighted density profile slope
$\gamma^{\prime}_{\mathrm{m}}$ (Equation 11) measured at
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ . Lower right: the relation between the
velocity dispersion profile outer slope and the density profile slope. In all
panels, The different ETGs with rising, flat, and dropping
$\gamma_{\mathrm{outer}}$ are indicated by the different-colored scattered
dots, following Fig. 4. The blue dotted-dashed lines and shaded regions stand
for the best linear fit and their [$16\%$, $84\%$] confidence intervals for
each of the correlations.

###  4.1 Correlation with density profile

Before proceeding into the details of exploring the density profiles, we make
a change in the quantity characterizing the non-Gaussian outer kinematic
structure from the kurtosis gradient $\Delta h_{4}(R)/\Delta\mathrm{log}R$ to
the outer kurtosis $\langle h_{4}\rangle$ (measured at
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}\pm 2.5\,\mathrm{kpc}$), which is the
luminosity-weighted $h_{4}$ for spaxels/Voronoi bins within a 5 kpc annulus
around $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. The motivation for our choice is
two fold: on the one hand, since we focus on the correlation with
$\gamma_{\mathrm{outer}}$, we would prefer to isolate the contribution from
outer $h_{4}$ alone; on the other hand, this change of variable can also avoid
the systematics in the stellar velocity distribution at small radii pertinent
to resolution effects (core softening) and feedback physics (dark matter
fraction) as discussed in Section 3.2. In the upper left panel of Fig. 7, we
show $\Delta h_{4}(R)/\Delta\mathrm{log}R$ versus the outer $h_{4}$ for the
IllustrisTNG ETGs. The best linear fit of the correlation has slope
$b=0.159\pm 0.022$, $p=5.12\times 10^{-12}$, and Pearson correlation
coefficient $r=0.456$, justifying that the outer $h_{4}$ is a decent
representative of the kurtosis gradient. Following that, we show in the upper
right panel of Fig. 7 that a positive correlation (slope $b=0.020\pm 0.012$,
$p=0.104$) also exists between the outer $h_{4}$ and the outer velocity
dispersion profile slope $\gamma_{\mathrm{outer}}$. Combining these two
panels, we anchor our following discussion regarding the link to density
profiles and assembly history on the outer $h_{4}$.

The power-law slope of the total density profile, $\gamma_{\rm PL}^{\prime}$
(same as $\gamma^{\prime}$ in Equation 1, renaming to avoid confusion with the
velocity dispersion profiles’ inner and outer slopes) reflects the underlying
gradient in circular velocity (interchangeable with velocity dispersion in the
isotropic case, Equation 2) such that $\gamma^{\prime}_{\mathrm{PL}}>2$
corresponds to a ‘dropping’ circular velocity profile,
$\gamma^{\prime}_{\mathrm{PL}}<2$ a ‘rising’ circular velocity profile, and a
flat circular velocity profile in the exact isothermal case
($\gamma^{\prime}_{\mathrm{PL}}=2$). In Paper I and II, this parameter was
obtained by a linear fit to the combined $\log\rho$-$\log r$ profile of dark
matter, stars, and gas in 100 logarithmic radial bins from
$0.4R_{\mathrm{eff}}$ to $4R_{\mathrm{eff}}$, where $R_{\mathrm{eff}}$ is the
projected stellar half mass radius. To study the total density profile at the
outskirts of these ETGs, we approximate the power-law slope with the ‘mass-
weighted’ logarithmic density slope (Dutton & Treu, 2014; Wang et al., 2020):

$\gamma^{\prime}_{\mathrm{m}}(R)=\frac{-1}{M(<R)}\int_{0}^{R}\gamma^{\prime}_{\mathrm{PL}}(r)4\pi
r^{2}\rho(r)dr=3-\frac{4\pi R^{3}\rho(R)}{M(<R)},$ (11)

where $M(<R)$ is the mass enclosed within a 3D spherical aperture with radius
$R$, and $\rho(R)$ is the average local density at radius $R$. We evaluate the
mass-weighted slope $\gamma^{\prime}_{\mathrm{m}}$ at
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ (local density $\rho(R)$ calculated in a
thin shell $\pm 2\%$ at this radius) to compare with the outer $h_{4}$ and
$\gamma_{\mathrm{outer}}$. The correlations along with the best linear fit
results are shown in the bottom panels in Fig. 7.

Figure 8: Left column: The relation between the outer $h_{4}$ values and
stellar assembly contributions. Major (minor) mergers are defined as mergers
having stellar mass ratio larger (smaller) than 0.25. Right column: The
relation between the outer $\beta_{\mathrm{3D}}$ and stellar assembly
contributions. In each column, the upper panel shows the link to the fraction
of stellar mass accreted in major mergers, while the lower panel shows the
link to the fraction of stellar mass accreted in minor perturbations, i.e.
minor mergers and flyby interactions. The key takeaway is that the outer
$\langle h_{4}\rangle$ correlates with minor perturbations but not major
mergers, while the velocity anisotropy $\langle\beta_{\mathrm{3D}}\rangle$
correlates more strongly with major mergers. The three types of ETGs with
rising, flat, and dropping outer velocity dispersion profiles are shown by the
scattered dots. The open markers with solid lines (dashed lines) denote the
median ($68\%$ scatter) of the rising, dropping, and flat
$\gamma_{\mathrm{outer}}$ ETGs combined. The blue dotted–dashed lines and
shaded regions stand for the best linear fit and their [$16\%$, $84\%$]
confidence intervals (please see text in Section 4.2 for best-fit slopes and
$p$-values.)

Intriguingly, both the outer $h_{4}$ (bottom left panel) and the velocity
dispersion profile outer slope $\gamma_{\mathrm{outer}}$ (bottom right panel)
are anti-correlated with the outer mass-weighted density profile slope at
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. Apparently, there is noticeable scatter
in both relations which lead to the scatter seen in the $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ correlation in Fig. 5.
Nevertheless, both of these negative correlations are statistically
significant, with $p=0.028$ for the outer
$h_{4}$-$\gamma^{\prime}_{\mathrm{m}}$ correlation (slope $b=-0.031\pm 0.014$)
and $p=0.027$ for the $\gamma_{\mathrm{outer}}$-$\gamma^{\prime}_{\mathrm{m}}$
correlation (slope $b=-0.177\pm 0.079$). With the velocity anisotropy ruled
out as a plausible explanation as discussed in Section 3.4, these two
correlations suggest that it is indeed the systematic variation of the total
density profile being shallower than isothermal
($\gamma^{\prime}_{\mathrm{m}}<2$) that leads to ‘rising’ gradients in the
velocity dispersion profiles ($\gamma_{\mathrm{outer}}>0$) and produces
positive outer $h_{4}$ in these ETGs (vice versa for the steeper than
isothermal case where $\gamma^{\prime}_{\mathrm{m}}>2$), giving rise to an
overall positive correlation between $\Delta h_{4}(R)/\Delta\mathrm{log}R$ and
$\gamma_{\mathrm{outer}}$. The physical interpretation behind this correlation
is that ongoing or past remnants of minor mergers producing positive outer
$h_{4}$ while simultaneously driving the total density profile shallower than
isothermal (Section 4.2). Moreover, combined with our finding in Paper I that
there is no significant correlation between the density profile slope and
velocity anisotropy (see Figure 8 in Wang et al. 2020), having the main
driving factor of a non-zero $h_{4}$ being a deviation from an isothermal
density profile allows for randomness between $h_{4}$ and
$\beta_{\mathrm{3D}}$, leading to the de-correlated
$h_{4}$-$\beta_{\mathrm{3D}}$ trend as seen in the bottom panel of Fig. 6.

###  4.2 Correlation with stellar assembly history

Galaxy mergers play an important role in the non-dissipative (dry) evolution
of ETGs at low redshift (Nipoti et al., 2009a; Nipoti et al., 2009b; Remus et
al., 2013; Wellons et al., 2016). In Paper I we demonstrated that the stellar
assembly history of ETGs parameterized by the fraction of in-situ-formed stars
($f_{\mathrm{in-situ}}$) positively correlates with the total density profile
power-law slope within $4R_{\mathrm{eff}}$ (see Figure 10 in Wang et al.
2020). In Paper II we have shown that gas-poor mergers at low redshift
($z\lesssim 0.5$), both major and minor mergers alike, work in coordination to
reduce $\gamma^{\prime}$ of the total density profile and establish the
positive $f_{\mathrm{in-situ}}-\gamma^{\prime}$ correlation (see Figure 7 in
Wang et al. 2019). Motivated by these finding and the positive relation of
$h_{4}$ at $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ and
$\gamma^{\prime}_{\mathrm{m}}$ as shown above in Fig. 7, we expect to also
find a correlation between the outer $h_{4}$ and the stellar assembly history
of these ETGs.

In Fig. 8, we present the relations between the outer $h_{4}$ (left column)
and outer $\beta_{\mathrm{3D}}$ (right column) measured at
$1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ as a function of the mass fraction of
stars in these ETGs accreted either from major mergers (upper row) or minor
perturbations combining the contribution from minor mergers and flyby
interactions (lower row). Major (minor) mergers are defined as mergers into
the main progenitors of our $z=0$ ETGs having stellar mass ratios larger
(smaller) than $1/4$. The solid line and open markers denote the median of the
combined sample including the rising, dropping, and flat outer velocity
dispersion profile ETGs, while the dashed lines represent the [$16\%$, $84\%$]
range ($1\,\sigma$) of the distributions. The outer $h_{4}$ shows almost no
correlation with the stellar constituent from major mergers (best linear fit
slope $b=-0.0006\pm 0.0103$, $p=0.9545$, statistically insignificant).
Intriguingly, the outer $h_{4}$ demonstrates a positive correlation (best
linear fit slope $b=0.0254\pm 0.0128$, $p=0.0474$, statistically significant)
with the stellar assembly contribution from minor perturbations at
$f_{\mathrm{minor\ merger\oplus flyby}}\lesssim 0.3$. This correlation is
shaped by the majority of galaxies in our sample which demonstrates consistent
increase in the median as well as the $1\,\sigma$ region. We also notice that
there is a slight decrease in outer $h_{4}$ at $f_{\mathrm{minor\ merger\oplus
flyby}}\sim 0.4$, beyond which $h_{4}$ remains largely constant. This
fluctuation in the outer $h_{4}$ at high minor perturbation fractions could be
due to the limited sample size at the massive end.

Nevertheless, the outer 3D velocity anisotropy shows positive correlations
with stellar assembly through both major mergers and minor perturbations, with
both relations having large scatter as well. The main difference regarding
these two correlations is that major mergers seem to drive a more significant
($\sim 0.4$ versus $\sim 0.1$) and steady (monotonic versus fluctuating)
increase in the outer $\beta_{\mathrm{3D}}$ compared to minor perturbations.
This is also evident from the steeper slope of the linear fit between
$\langle\beta_{\mathrm{3D}}\rangle$ and $f_{\mathrm{major\ merger}}$ (slope
$b=0.702\pm 0.115$, $p=5.83\times 10^{-9}$, very significant) versus
$f_{\mathrm{minor\ merger\oplus flyby}}$ (slope $b=0.196\pm 0.157$, $p=0.213$,
poor statistical significance). Hence, the positive outer $h_{4}$ in ETGs
tends to be built up by minor mergers and galaxy interactions (at least for
the majority of galaxies with $f_{\mathrm{minor\ merger\oplus flyby}}\lesssim
0.4$), but the radial anisotropy of stars ($\beta_{\mathrm{3D}}$) at the
galaxy outskirts is mainly driven by major mergers, resulting in an almost
decoupled $h_{4}$-$\beta_{\mathrm{3D}}$ relation as seen in Fig. 6.

Although the four sets of relations shown above all possess scatter, their
differences hint that the seemingly decoupled outer $h_{4}$ and
$\beta_{\mathrm{3D}}$ (Fig. 6) trends originate from the different impacts of
major and minor perturbations during the formation process of ETGs. Major
mergers can create cuspy (wet merger) or cored (dry merger) central stellar
densities depending on the gas fraction of the merger (Hopkins et al., 2008c,
2009a; Kormendy et al., 2009; Hopkins et al., 2009b, d; Tacchella et al.,
2019), modifying the central density profile. Major mergers can also
drastically change the morphology and dynamics of a massive galaxy from being
tangentially biased disk galaxies (late-type) to radial-orbit-dominated ETGs
(Kormendy et al., 2009; Kormendy & Bender, 2012; Rodriguez-Gomez et al.,
2017), which explains the positive correlation between $\beta_{\mathrm{3D}}$
and $f_{\mathrm{major-merger}}$. Although in isolated galaxy merger
simulations, a bump in the velocity dispersion profile can also be created at
radii $>R_{\mathrm{eff}}$ via major mergers (Schauer et al., 2014),
simultaneous minor mergers could counter-act that effect in a cosmological
setting ($1\leqslant z\leqslant 2$, Wang et al. 2019). In addition, minor
mergers contribute significantly to the growth of ETGs compared to major
mergers (Guo & White, 2008), perpetuate the inside-out growth of extensive
stellar halos and modify stellar velocity distributions of ETGs at large radii
(Murante et al., 2007; Bernardi, 2009; Khochfar et al., 2011; Hilz et al.,
2012; Hilz et al., 2013; Pulsoni et al., 2021; Dolfi et al., 2021).

By visually inspecting the kinematic maps of $h_{3}$ and $h_{4}$ for
individual galaxies, we find that at the locations where $h_{3}$ and $h_{4}$
locally peak, it is also a common location to find a coinciding overdensity in
the projected luminosity map (see Appendix B for examples). These
overdensities reflect the transient infall phases of accreted satellite
galaxies on to the host ETGs, before they are tidally-disrupted and identified
as minor mergers by Sublink. The corresponding LOS velocity profile in the
spaxels/Voronoi-bins that cut through these satellites show peaks that are
phase from the mean velocity of the galaxy, which leads to large $h_{3}$ and
$h_{4}$ values in those spaxels/Voronoi-bins. Eventually, the cumulative
effect of these minor perturbations (which are predominantly dry mergers at
low redshift, see Figure 7 Wang et al. 2019 for simulation evidence and
Derkenne et al. 2021 for latest dynamical modeling constraints from VLT-MUSE)
not only leads to positive outer $h_{4}$, but also drives the outer density
profile to be shallower than isothermal ($\gamma^{\prime}_{\mathrm{m}}$),
establishing the positive correlation between $\Delta
h_{4}(R)/\Delta\mathrm{log}R$ and $\gamma_{\mathrm{outer}}$. Combined with the
aforementioned impact of minor mergers on the formation of ETGs, the formation
of the non-Gaussian outer kinematic structure in early-type galaxies is a
natural consequence of hierarchical structure formation, and it co-evolves
with the outer density profile mainly via contributions from minor
perturbations. As we will discuss in the following, these trends are also
consistent with the results from the probes of galaxy (halo) environment.

###  4.3 Correlation with environment

Figure 9: Left panel: The correlation between the outer velocity dispersion
profile slope ($\gamma_{\mathrm{outer}}$) and the total stellar mass
($M_{\ast}$) of our simulated ETGs, color-coded by their outer $h_{4}$ values
measured at around $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. Right panel: The
correlation between $\gamma_{\mathrm{outer}}$ and the host halo mass
($M_{200}$) of our simulated ETGs, color-coded by their outer $h_{4}$ values
measured at around $1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$. The open markers with
solid lines in both panels denote the median, while the dashed lines represent
the [$16\%$, $84\%$] of the distribution for $\gamma_{\mathrm{outer}}$ in
different stellar mass and halo mass bins. The red crossed with dotted-dashed
lines are median values for the same correlations from Veale et al. (2018).
The upper side plots show the fraction of galaxies with rising (blue), flat
(grey), and dropping (red) outer dispersion profiles in each stellar/halo mass
bin. The error bars in each bin stand for the $68\%$ Wilson score confidence
intervals.

Another interesting feature related to the $\Delta
h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$ found in Veale et al.
(2018) is that there is an increasing fraction of ETGs with ‘rising’ outer
velocity dispersion profiles living in denser environments and more massive
halos. Again, the authors speculate that this is also induced by the
systematic variation of density profiles instead of velocity anisotropy across
different ETGs. Since the environmental measures $\nu_{10}$ and $\delta_{g}$
in Veale et al. (2017); Veale et al. (2018) are based on the galaxy K-band
luminosity weighted on the MASSIVE Survey limiting magnitude and completeness,
we do not attempt to perform mock observations of our simulated ETGs based on
raw SDSS $K$-band magnitude output by IllustrisTNG. Hence, we focus on
comparing the total stellar ($M_{\ast}$) and host dark matter halo masses
($M_{200}$, mass contained in the radius where the mean density is $200\times$
the universal critical density) of our ETG sample to the trends observed in
the MASSIVE Survey. The stellar and halo masses serve as sensible proxies of
the ETG environment as all targets in MASSIVE are central galaxies and we
follow the same selection for our mock IllustrisTNG ETGs.

In Fig. 9, we show $\gamma_{\mathrm{outer}}$, as well as the fraction of ETGs
having ‘rising’, ‘flat’, or ‘dropping’ velocity dispersion profile outer
slopes, against their stellar (left panel) and halo masses (right panel). The
scattered dots are individual galaxies colored by their outer $h_{4}$ values,
the solid curve with diamond markers denotes the binned median, and the dashed
lines represent the [$16\%$, $84\%$] interval of the distribution within each
mass bin. We see that $\gamma_{\mathrm{outer}}$ increases with both increasing
stellar and halo masses, in agreement with the trend of MASSIVE median values
(red crosses). We notice that the MASSIVE ETGs (89 galaxies, Veale et al.
2018, also see their Figures 7, 8, and 9 which also shows large scatter in
$\gamma_{\mathrm{outer}}$) have a slightly steeper relation with $M_{\ast}$,
possibly due to a more complete sample of heavy galaxies with
$\log_{10}(M_{\ast}/\mathrm{M_{\astrosun}})\geqslant 11.6$ compared to the
IllustrisTNG ETG sample (18 galaxies). The same trend is reflected in the
upper side plots for the fraction of galaxies with different velocity
dispersion profile outer slopes. The general trend in both columns indicate
that higher fractions of galaxies will have ‘rising’ $\gamma_{\mathrm{outer}}$
with increasing galaxy and halo mass (the error bars in the upper panels are
the Wilson score $68\%$ confidence intervals, Wilson 1927).

Recasting on the findings in Paper I (Wang et al., 2020), the power-law slope
of the total density profile anti-correlating with the total stellar mass
(Figure 4) is driven by the variation in the inner dark matter density profile
(Figure 17), and hence anti-correlating with the halo mass (Figure 15). Given
that the variation in $\gamma_{\mathrm{outer}}$ is mainly driven by a negative
correlation with the density profile (lower right panel in Fig. 7), the
findings in Paper I self-consistently predict the trends seen in Fig. 9, once
again confirming the conjecture for the observed phenomena in Veale et al.
(2018). The color index in the figure shows that the outer $h_{4}$ also tend
to be higher in ETGs with higher stellar and halo mass, consistent with the
positive $\Delta h_{4}(R)/\Delta\mathrm{log}R$-$\gamma_{\mathrm{outer}}$
relation.

The physical interpretation for the emergence of the
$\gamma_{\mathrm{outer}}$-mass correlation also fits in to the multi-phase
formation path of near-isothermal ETG density profiles in Paper II (Wang et
al., 2019): heavier ETGs are quenched earlier by AGN feedback (Figures 4 and
8) and enter earlier into the non-dissipative evolution phase dominated by
gas-poor mergers at lower redshift $z\lesssim 2$ (Figure 7), where the total
density profile progressively evolves shallower, and more massive galaxies
living in more massive halos (denser environments) experience more mergers.
Although ETGs go through both major and minor mergers, with previous studies
suggesting that major mergers start to dominate the ex-situ stellar population
at stellar masses $\log_{10}(M_{\ast}/\mathrm{M_{\astrosun}})\gtrsim 11.2$
(e.g., Fig. 10 in Tacchella et al. 2019), which is also the case for our ETG
sample with median $f_{\mathrm{major\ merger}}$ = 0.31 and median
$f_{\mathrm{minor\ merger\oplus flyby}}=0.12$, the stellar mass fraction from
both minor and major mergers increase with $M_{\ast}$. Since the outer $h_{4}$
is mainly driven by minor perturbations instead of major mergers (Section 4.2,
lower left panel in Fig. 8), the mass trend of $f_{\mathrm{minor\ merger\oplus
flyby}}$ alone dictates the mass trend of $\gamma_{\mathrm{outer}}$ and outer
$h_{4}$. This eventually creates the positive correlation between
$\gamma_{\mathrm{outer}}$ and the kurtosis gradient through the mass-dependent
variation of the total density profile, indicating that the formation of the
outer (non-Gaussian) kinematic structure in ETGs stands as a natural
consequence of hierarchical galaxy-halo assembly.

##  5 Conclusions and Outlook

We study the relation between the kinematic structure and density profiles of
207 massive
($\log_{10}\left(M_{\ast}/\mathrm{M_{\astrosun}}\right)\in[11.0,12.4]$, non-
major-merger, Fig. 2) early-type-galaxies selected from the IllustrisTNG
Simulations (Marinacci et al., 2018; Naiman et al., 2018; Nelson et al., 2018;
Pillepich et al., 2018b; Springel et al., 2018). We produce mock integral-
field-unit spectroscopic observations (Fig. 1) for these ETGs in random
projections and compute the second (velocity dispersion $\sigma$), third
(skewness $h_{3}$), and fourth (kurtosis $h_{4}$) velocity moments of the
projected stellar velocity profiles (Fig.3). The outer slope
($\gamma_{\mathrm{outer}}$ ) of the LOSVD profiles (Fig. 4) and the key
positive correlation between $\gamma_{\mathrm{outer}}$ and the kurtosis
gradient ($\Delta h_{4}(R)/\Delta\mathrm{log}R$, Fig. 5) for the simulated
ETGs are consistent with that found in the MASSIVE Survey (Ma et al., 2014;
Veale et al., 2017; Veale et al., 2018). This realistic mock ETG sample serves
as a benchmark for disentangling the correlation between the outer non-
Gaussian kinematic features and velocity dispersion profiles in observed ETGs
from the MASSIVE Survey (Veale et al., 2018). In particular, we find:

  * •

The average velocity anisotropy and $h_{4}$ show almost _no_ correlation with
the 3D velocity anisotropy $\beta_{\mathrm{3D}}$ at radii $\geqslant
R_{\mathrm{eff}}^{\mathrm{MGE}}$ (Fig. 6 lower panel), indicating that
positive outer $h_{4}$ values is unlikely to be produced by stellar velocity
anisotropy under fixed density profiles (Gerhard, 1993).

  * •

The outer $h_{4}$ measured at $\sim 1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ is a
good proxy for the kurtosis gradient and also positively correlates with
$\gamma_{\mathrm{outer}}$ (Fig. 7 upper row). The mass-weighted slope of the
total density profile at that radius negatively correlates with both the outer
$h_{4}$ and $\gamma_{\mathrm{outer}}$ (Fig. 7 lower row). This justifies the
conjecture (Veale et al., 2018) that systematic variations around an exact
isothermal profile drive the positive correlation between $\Delta
h_{4}(R)/\Delta\mathrm{log}R$ and $\gamma_{\mathrm{outer}}$.

  * •

The outer $h_{4}$ positively correlates with the mass fraction of stars
accreted via minor perturbations (minor mergers and flyby interactions), while
the outer velocity anisotropy correlates more significantly with the stellar
mass fraction from major mergers (Fig. 8), resulting in a seemingly decoupled
$h_{4}$–$\beta_{\mathrm{3D}}$ relation at the ETG outskirts (Fig. 6).

  * •

The values of $\gamma_{\mathrm{outer}}$, outer $h_{4}$, and the fraction of
ETGs with ‘rising’ outer velocity dispersion profiles increase with increasing
stellar and halo masses. This is mainly driven by the increase in minor merger
fractions with increasing ETG stellar mass affecting the outer kinematic
structure and density profile simultaneously (Fig. 9).

As an augmentation to the detailed studies of the ETG density profiles in
Paper I and II, this work generalizes those findings and links them with the
outer kinematic structure in these galaxies. Our finding that the outer
$h_{4}$ correlates with minor mergers and the density profile, while the outer
$\beta_{\mathrm{3D}}$ correlates with major mergers, also explains the
apparent stochasticity between the total density profile and stellar velocity
anisotropy of ETGs found in Figure 8 of Paper I. The outer $h_{4}$–minor
merger correlation is also consistent with our previous finding: in systems
where the stellar assembly in ETGs is dominated by minor perturbations at low
redshift ($z<0.5$), a higher accreted stellar fraction leads to steeper
density profile slopes (Figure 10, Paper I). The dependence of
$\gamma_{\mathrm{outer}}$ and outer $h_{4}$ with environment is consistent
with the previous findings that both stellar and halo masses anti-correlate
with the steepness of the total density profile (Figures 4, 15, and 17 in
Paper I), and fits in smoothly to the low-redshift ($z<0.5$) minor-merger-
driven evolution phase of the ETG density profiles as discussed in Paper II
(Figures 4, 7, 8). Combining the effects of stellar assembly on outer $h_{4}$,
we conclude that the outer kinematic structure of ETGs co-evolves with their
total density profiles, which is a natural consequence of hierarchical
assembly in cosmological structure formation. The broad agreement of our
results with the MASSIVE Survey ETGs also highlights the power of IllustrisTNG
to elucidate the underlying correlations in a realistic formation scenario,
that could not otherwise be discerned from observations alone.

Apart from the many aspects of our findings that do reproduce and explain the
observed ETG kinematic features, we highlight that there are important
systematics involving the central stellar velocity dispersion of the
IllustrisTNG galaxies. We observe more galaxies with positive
$\gamma_{\mathrm{inner}}$ compared to MASSIVE ETGs (Fig. 4) indicating
underestimation of the central stellar velocity dispersion, which is
consistent with the tendency of over-predicting central dark matter fraction
(producing halo contraction) in IllustrisTNG galaxies (Lovell et al., 2018;
Wang et al., 2020). This is also related to the negative
$h_{4}$–$\beta_{\mathrm{3D}}$ correlation within
$0.5R_{\mathrm{eff}}^{\mathrm{MGE}}$ as seen in Fig. 6. These limitations pave
a way forward for better understanding the interplay between dark matter and
baryons, as well as for future improvements of the simulation subgrid physics
model (especially related to AGN feedback), which could be better constrained
by refined kinematic features at both large and small galactic radii.

Another interesting feature is that the presence of infalling substructures
can significantly boost $h_{3}$ and $h_{4}$ along the line-of-sight at the
outskirts of ETGs (Appendix B). Although quantifying the impact of mergers and
their remnants on the kinematic structure of ETGs is beyond the scope of this
work, this finding ties in neatly to our finding of the ETG outer kinematic
structure by the fraction of stars acquired through minor perturbations. With
ultra-high resolution IFU spectroscopy from MUSE (Bacon et al., 2010) on the
Very Large Telescope (VLT) or KCWI (Morrissey et al., 2018) on the Keck
Telescope, non-Gaussian kinematic features could potentially provide a novel
approach for discovering faint satellites around ETG hosts. Combined with
upcoming photometric campaigns such as the Vera C. Rubin Observatory (LSST,
Ivezić et al. 2019), the Nancy Grace Roman Space Telescope (WFIRST, Spergel et
al. 2015), and Euclid (Amiaux et al., 2012) that advance in survey depth and
field of view, we will soon be able to establish a more profound understanding
of the co-evolution between the kinematic structure, density profile, and even
the build up of the stellar halo around ETGs in the broad context of
hierarchical structure formation.

## Acknowledgements

We thank Dandan Xu, Ethan O. Nadler, Philip Mansfield, and Shengdong Lu for
helpful discussions and support during the preparation of this paper. We thank
the anonymous referee for insightful comments that helped to improve the
draft. YW acknowledges the past support from the Tsinghua Xuetang Talents
Programme and current support from a Stanford-KIPAC Chabolla Fellowship. This
work is partly supported by the National Key Basic Research and Development
Programme of China (No. 2018YFA0404501 to SM), and by the National Science
Foundation of China (Grant No. 11821303 and 11761131004 to SM). MV
acknowledges support through an MIT RSC award, a Kavli Research Investment
Fund, NASA ATP grant NNX17AG29G, and NSF grants AST-1814053 and AST-1814259.
VS acknowledges support by the Deutsche Forschungsgemeinschaft through project
SP 709/5-1. This research made use of computational resources at the
MIT/Harvard computing facilities supported by FAS and MIT MKI; the authors are
thankful for the support from the FAS Research Computing team. The flagship
simulations of the IllustrisTNG project used in this work have been run on the
HazelHen Cray XC40-system at the High Performance Computing Center Stuttgart
as part of project GCS-ILLU of the Gauss Centre for Supercomputing (GCS). This
research made extensive use of arXiv.org and NASA’s Astrophysics Data System
for bibliographic information.

## Data Availability

The IllustrisTNG public data can be accessed at https://www.tng-
project.org/data/. The data in this article that are not part of the
IllustrisTNG data release will be shared upon reasonable request to the
corresponding author.

## References

  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
  * Amiaux et al. (2012) Amiaux J., et al., 2012, in Clampin M. C., Fazio G. G., MacEwen H. A., Oschmann Jacobus M. J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave. p. 84420Z (arXiv:1209.2228), doi:10.1117/12.926513
  * Auger et al. (2010) Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Moustakas L. A., Burles S., 2010, ApJ, 724, 511
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508, doi:10.1117/12.856027
  * Baes et al. (2005) Baes M., Dejonghe H., Buyle P., 2005, A&A, 432, 411
  * Barnabè et al. (2009) Barnabè M., Czoske O., Koopmans L. V. E., Treu T., Bolton A. S., Gavazzi R., 2009, MNRAS, 399, 21
  * Barnabè et al. (2011) Barnabè M., Czoske O., Koopmans L. V. E., Treu T., Bolton A. S., 2011, MNRAS, 415, 2215
  * Bell & de Jong (2001) Bell E. F., de Jong R. S., 2001, ApJ, 550, 212
  * Bellstedt et al. (2018) Bellstedt S., et al., 2018, MNRAS, 476, 4543
  * Bender et al. (1992) Bender R., Burstein D., Faber S. M., 1992, ApJ, 399, 462
  * Bernardi (2009) Bernardi M., 2009, MNRAS, 395, 1491
  * Binney & Mamon (1982) Binney J., Mamon G. A., 1982, MNRAS, 200, 361
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Cappellari (2002) Cappellari M., 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari M., 2008, MNRAS, 390, 71
  * Cappellari et al. (2015) Cappellari M., et al., 2015, ApJ, 804, L21
  * Cassata et al. (2013) Cassata P., et al., 2013, ApJ, 775, 106
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Cole et al. (2000) Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
  * De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
  * De Lucia et al. (2006) De Lucia G., Springel V., White S. D. M., Croton D., Kauffmann G., 2006, MNRAS, 366, 499
  * Dejonghe & Merritt (1992) Dejonghe H., Merritt D., 1992, ApJ, 391, 531
  * Derkenne et al. (2021) Derkenne C., McDermid R. M., Poci A., Remus R.-S., Jørgensen I., Emsellem E., 2021, MNRAS, 506, 3691
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  * Diemer et al. (2019) Diemer B., et al., 2019, MNRAS, 487, 1529
  * Djorgovski & Davis (1987) Djorgovski S., Davis M., 1987, ApJ, 313, 59
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Dolfi et al. (2021) Dolfi A., Forbes D. A., Couch W. J., Bekki K., Ferré-Mateu A., Romanowsky A. J., Brodie J. P., 2021, MNRAS, 504, 4923
  * Donnari et al. (2019) Donnari M., et al., 2019, MNRAS, 485, 4817
  * Dressler (1980) Dressler A., 1980, ApJ, 236, 351
  * Dutton & Treu (2014) Dutton A. A., Treu T., 2014, MNRAS, 438, 3594
  * Emsellem et al. (1994) Emsellem E., Monnet G., Bacon R., Nieto J. L., 1994, A&A, 285, 739
  * Ene et al. (2019) Ene I., Ma C.-P., McConnell N. J., Walsh J. L., Kempski P., Greene J. E., Thomas J., Blakeslee J. P., 2019, ApJ, 878, 57
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Gavazzi et al. (2007) Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, ApJ, 667, 176
  * Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  * Genel et al. (2018) Genel S., et al., 2018, MNRAS, 474, 3976
  * Gerhard (1993) Gerhard O. E., 1993, MNRAS, 265, 213
  * Gerhard et al. (1998) Gerhard O., Jeske G., Saglia R. P., Bender R., 1998, MNRAS, 295, 197
  * Greene et al. (2019) Greene J. E., et al., 2019, ApJ, 874, 66
  * Guo & White (2008) Guo Q., White S. D. M., 2008, MNRAS, 384, 2
  * Habouzit et al. (2019) Habouzit M., et al., 2019, MNRAS, 484, 4413
  * Hemler et al. (2021) Hemler Z. S., et al., 2021, MNRAS, 506, 3024
  * Hill et al. (2008) Hill G. J., MacQueen P. J., Palunas P., Barnes S. I., Shetrone M. D., 2008, in McLean I. S., Casali M. M., eds, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, pp 92 – 106, doi:10.1117/12.788225, https://doi.org/10.1117/12.788225
  * Hilz et al. (2012) Hilz M., Naab T., Ostriker J. P., Thomas J., Burkert A., Jesseit R., 2012, MNRAS, 425, 3119
  * Hilz et al. (2013) Hilz M., Naab T., Ostriker J. P., 2013, MNRAS, 429, 2924
  * Hopkins et al. (2006) Hopkins P. F., Hernquist L., Cox T. J., Robertson B., Springel V., 2006, ApJS, 163, 50
  * Hopkins et al. (2008a) Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008a, ApJS, 175, 356
  * Hopkins et al. (2008b) Hopkins P. F., Cox T. J., Kereš D., Hernquist L., 2008b, ApJS, 175, 390
  * Hopkins et al. (2008c) Hopkins P. F., Hernquist L., Cox T. J., Dutta S. N., Rothberg B., 2008c, ApJ, 679, 156
  * Hopkins et al. (2009a) Hopkins P. F., Cox T. J., Dutta S. N., Hernquist L., Kormendy J., Lauer T. R., 2009a, ApJS, 181, 135
  * Hopkins et al. (2009b) Hopkins P. F., Lauer T. R., Cox T. J., Hernquist L., Kormendy J., 2009b, ApJS, 181, 486
  * Hopkins et al. (2009c) Hopkins P. F., et al., 2009c, MNRAS, 397, 802
  * Hopkins et al. (2009d) Hopkins P. F., Cox T. J., Younger J. D., Hernquist L., 2009d, ApJ, 691, 1168
  * Huertas-Company et al. (2013) Huertas-Company M., et al., 2013, MNRAS, 428, 1715
  * Humphrey & Buote (2010) Humphrey P. J., Buote D. A., 2010, MNRAS, 403, 2143
  * Humphrey et al. (2006) Humphrey P. J., Buote D. A., Gastaldello F., Zappacosta L., Bullock J. S., Brighenti F., Mathews W. G., 2006, ApJ, 646, 899
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Johansson et al. (2012) Johansson P. H., Naab T., Ostriker J. P., 2012, ApJ, 754, 115
  * Khochfar et al. (2011) Khochfar S., et al., 2011, MNRAS, 417, 845
  * King (2003) King A., 2003, ApJ, 596, L27
  * Koopmans et al. (2006) Koopmans L. V. E., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, ApJ, 649, 599
  * Koopmans et al. (2009) Koopmans L. V. E., et al., 2009, ApJ, 703, L51
  * Kormendy & Bender (2012) Kormendy J., Bender R., 2012, ApJS, 198, 2
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kormendy et al. (2009) Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216
  * Krajnović et al. (2011) Krajnović D., et al., 2011, MNRAS, 414, 2923
  * Li et al. (2016) Li H., Li R., Mao S., Xu D., Long R. J., Emsellem E., 2016, MNRAS, 455, 3680
  * Li et al. (2018a) Li H., et al., 2018a, MNRAS, 476, 1765
  * Li et al. (2018b) Li R., Shu Y., Wang J., 2018b, MNRAS, 480, 431
  * Li et al. (2019) Li R., et al., 2019, MNRAS, 490, 2124
  * Lovell et al. (2018) Lovell M. R., et al., 2018, MNRAS, 481, 1950
  * Lyskova et al. (2018) Lyskova N., Churazov E., Naab T., 2018, MNRAS, 475, 2403
  * Ma et al. (2014) Ma C.-P., Greene J. E., McConnell N., Janish R., Blakeslee J. P., Thomas J., Murphy J. D., 2014, ApJ, 795, 158
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * McGaugh et al. (2000) McGaugh S. S., Schombert J. M., Bothun G. D., de Blok W. J. G., 2000, ApJ, 533, L99
  * Merritt & Saha (1993) Merritt D., Saha P., 1993, ApJ, 409, 75
  * Morrissey et al. (2018) Morrissey P., et al., 2018, ApJ, 864, 93
  * Moster et al. (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Murante et al. (2007) Murante G., Giovalli M., Gerhard O., Arnaboldi M., Borgani S., Dolag K., 2007, MNRAS, 377, 2
  * Naab et al. (2007) Naab T., Johansson P. H., Ostriker J. P., Efstathiou G., 2007, ApJ, 658, 710
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelson et al. (2015) Nelson D., et al., 2015, Astronomy and Computing, 13, 12
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Nipoti et al. (2009a) Nipoti C., Treu T., Bolton A. S., 2009a, ApJ, 703, 1531
  * Nipoti et al. (2009b) Nipoti C., Treu T., Auger M. W., Bolton A. S., 2009b, ApJ, 706, L86
  * Nipoti et al. (2020) Nipoti C., Cannarozzo C., Calura F., Sonnenfeld A., Treu T., 2020, MNRAS, 499, 559
  * Oser et al. (2010) Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312
  * Paturel et al. (2003) Paturel G., Petit C., Prugniel P., Theureau G., Rousseau J., Brouty M., Dubois P., Cambrésy L., 2003, A&A, 412, 45
  * Pillepich et al. (2018a) Pillepich A., et al., 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich A., et al., 2018b, MNRAS, 475, 648
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Poci et al. (2017) Poci A., Cappellari M., McDermid R. M., 2017, MNRAS, 467, 1397
  * Pulsoni et al. (2021) Pulsoni C., Gerhard O., Arnaboldi M., Pillepich A., Rodriguez-Gomez V., Nelson D., Hernquist L., Springel V., 2021, A&A, 647, A95
  * Read & Steger (2017) Read J. I., Steger P., 2017, MNRAS, 471, 4541
  * Remus et al. (2013) Remus R.-S., Burkert A., Dolag K., Johansson P. H., Naab T., Oser L., Thomas J., 2013, ApJ, 766, 71
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez V., et al., 2016, MNRAS, 458, 2371
  * Rodriguez-Gomez et al. (2017) Rodriguez-Gomez V., et al., 2017, MNRAS, 467, 3083
  * Romanowsky et al. (2003) Romanowsky A. J., Douglas N. G., Arnaboldi M., Kuijken K., Merrifield M. R., Napolitano N. R., Capaccioli M., Freeman K. C., 2003, Science, 301, 1696
  * Ruff et al. (2011) Ruff A. J., Gavazzi R., Marshall P. J., Treu T., Auger M. W., Brault F., 2011, ApJ, 727, 96
  * Schauer et al. (2014) Schauer A. T. P., Remus R.-S., Burkert A., Johansson P. H., 2014, ApJ, 783, L32
  * Serra et al. (2016) Serra P., Oosterloo T., Cappellari M., den Heijer M., Józsa G. I. G., 2016, MNRAS, 460, 1382
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Silk & Rees (1998) Silk J., Rees M. J., 1998, A&A, 331, L1
  * Sonnenfeld et al. (2013) Sonnenfeld A., Treu T., Gavazzi R., Suyu S. H., Marshall P. J., Auger M. W., Nipoti C., 2013, ApJ, 777, 98
  * Sonnenfeld et al. (2014) Sonnenfeld A., Nipoti C., Treu T., 2014, ApJ, 786, 89
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stevens et al. (2019) Stevens A. R. H., et al., 2019, MNRAS, 483, 5334
  * Tacchella et al. (2019) Tacchella S., et al., 2019, MNRAS, 487, 5416
  * Terrazas et al. (2020) Terrazas B. A., et al., 2020, MNRAS, 493, 1888
  * Torrey et al. (2014) Torrey P., Vogelsberger M., Genel S., Sijacki D., Springel V., Hernquist L., 2014, MNRAS, 438, 1985
  * Torrey et al. (2018) Torrey P., et al., 2018, MNRAS, 477, L16
  * Torrey et al. (2019) Torrey P., et al., 2019, MNRAS, 484, 5587
  * Tortora et al. (2014) Tortora C., La Barbera F., Napolitano N. R., Romanowsky A. J., Ferreras I., de Carvalho R. R., 2014, MNRAS, 445, 115
  * Veale et al. (2017) Veale M., et al., 2017, MNRAS, 464, 356
  * Veale et al. (2018) Veale M., Ma C.-P., Greene J. E., Thomas J., Blakeslee J. P., Walsh J. L., Ito J., 2018, MNRAS, 473, 5446
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
  * Vogelsberger et al. (2014a) Vogelsberger M., et al., 2014a, MNRAS, 444, 1518
  * Vogelsberger et al. (2014b) Vogelsberger M., et al., 2014b, Nature, 509, 177
  * Vogelsberger et al. (2018) Vogelsberger M., et al., 2018, MNRAS, 474, 2073
  * Vogelsberger et al. (2020a) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020a, Nature Reviews Physics, 2, 42
  * Vogelsberger et al. (2020b) Vogelsberger M., et al., 2020b, MNRAS, 492, 5167
  * Wang et al. (2019) Wang Y., et al., 2019, MNRAS, 490, 5722
  * Wang et al. (2020) Wang Y., et al., 2020, MNRAS, 491, 5188
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Wellons et al. (2015) Wellons S., et al., 2015, MNRAS, 449, 361
  * Wellons et al. (2016) Wellons S., et al., 2016, MNRAS, 456, 1030
  * Williams et al. (2010) Williams M. J., Bureau M., Cappellari M., 2010, MNRAS, 409, 1330
  * Wilson (1927) Wilson E. B., 1927, Journal of the American Statistical Association, 22, 209
  * Wyithe & Loeb (2003) Wyithe J. S. B., Loeb A., 2003, ApJ, 595, 614
  * Xie et al. (2015) Xie L., Guo Q., Cooper A. P., Frenk C. S., Li R., Gao L., 2015, MNRAS, 447, 636
  * Xu et al. (2017) Xu D., Springel V., Sluse D., Schneider P., Sonnenfeld A., Nelson D., Vogelsberger M., Hernquist L., 2017, MNRAS, 469, 1824
  * Xu et al. (2019) Xu D., et al., 2019, MNRAS, 489, 842
  * Yun et al. (2019) Yun K., et al., 2019, MNRAS, 483, 1042
  * van der Marel & Franx (1993) van der Marel R. P., Franx M., 1993, ApJ, 407, 525

##  Appendix A The effective radius

Figure 10: The comparison between the projected 2D effective radius calculated
with all stellar particles assigned to the galaxy determined by Subfind,
versus the MGE effective radius output by our mock observation pipeline. The
MGE effective radii closely follow the 2D projected effective radii for
galaxies with small sizes ($R_{\mathrm{eff}}\lesssim 10$ kpc). For the larger
galaxies, the MGE effective radii are significantly smaller than the projected
effective radii, which results from the MGE approach removing most of the
diffuse intra-cluster light component assigned to the galaxy by Subfind,
giving more realistic estimates of the effective radii for larger galaxies.
Figure 11: Three examples of spaxels/Voronoi bins cutting through infalling
satellite galaxies around two different host galaxies that introduces large
$h_{3}$ and $h_{4}$. Upper panel: Host Subfind ID 69507, with stellar mass
$10^{12.35}\mathrm{M_{\astrosun}}$ and halo mass
$10^{14.32}\mathrm{M_{\astrosun}}$. Lower panel: Host Subfind ID 257302, with
stellar mass $10^{11.59}\mathrm{M_{\astrosun}}$ and halo mass
$10^{13.32}\mathrm{M_{\astrosun}}$. The dashed circles in the luminosity maps
indicate the location of these satellites, while the red crosses indicate the
location of the spaxels/Voronoi bins covering them. Blue histograms show the
LOS stellar velocity distribution in those bins (centered on the mean
velocity, $v=0$, large peak offsets from the zero point indicate satellite
stars), whereas the red curves indicate the $4^{\mathrm{th}}$ order Gauss-
Hermite function fit. Their corresponding locations in the $h_{3}$ and $h_{4}$
maps are also marked by black circles.

In this section, we show the differences between the projected 2D effective
radius ($R_{\mathrm{eff}}^{\mathrm{2D,proj}}$) and
$R_{\mathrm{eff}}^{\mathrm{MGE}}$ of galaxies given by the best MGE fit.
$R_{\mathrm{eff}}^{\mathrm{2D,proj}}$ is measured by projecting galaxies in
random directions (along the z axis of the simulation box) and searching for
the radius at which the enclosed projected stellar mass is half of the total
stellar mass. However, this procedure includes all stellar particles assigned
to galaxies by Subfind, and includes large contributions from the intra-
cluster light (ICL) component around the massive ETGs that we selected.
Another approach is to first model the galaxy’s projected luminosity
distribution in its central region (square aperture with 80 kpc side length)
using Multi-Gaussian Expansion (Emsellem et al., 1994), and calculate the
half-light radius using the set of analytic best-fit 2D Gaussians which
eventually removes most of the ICL. It also has the advantage over simpler
effective radius definitions in 2MASS (semi major axis length of an ellipse
enclosing half the total galaxy luminosity) and the NSA catalog (half light
radius from a single S$\mathrm{\acute{e}}$rsic fit) to flexibly preserve
information of the spatial distribution for different components in the galaxy
luminosity profile.

The comparison between the effective radii derived using these two methods is
shown in Fig. 10. For smaller galaxies ($R_{\mathrm{eff}}\lesssim 10$ kpc),
the two definitions yield consistent size measurements, while for larger
galaxies the MGE effective radii are significantly smaller than the 2D
projected effective radii, which is due to the removal of the ICL. We adopt
$R_{\mathrm{eff}}^{\mathrm{MGE}}$ in our comparison with observational data.

##  Appendix B Examples of minor mergers impacting the outer kinematic
structure

In this section we showcase three examples of ongoing minor mergers
introducing large non-Gaussian moments to the outer kinematic structure of
ETGs. In Fig. 11, we mark out three infalling satellites (overdensities in the
luminosity maps) around two different host ETGs, with spaxels/Voronoi bins
cutting through them along the line-of-sight, and resulting in significant
$h_{3}$ and $h_{4}$ values. Clearly, these three pixels all demonstrate large
peak offsets from the mean velocity (centered on $v=0$). While $h_{3}$
produced in this manner may diminish over time as the satellite orbits around
the host due to its odd parity, a positive $h_{4}$ can be preserved over time
due to its even parity. The cumulative effect of such in many minor mergers
(along with flyby interactions that tidally heat galaxies, accelerating stars
into the tails of the LOS velocity distribution) can gradually build up a
global positive $h_{4}$ in the outskirts ($\gtrsim
1.5R_{\mathrm{eff}}^{\mathrm{MGE}}$) of ETGs, leading to the positive
correlation of the outer $h_{4}$ and stellar mass fraction from minor
perturbations (minor mergers plus flyby interactions) as seen in Fig. 7.